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Abstract 

The basic objective of this research is to extend the capabilities of Large Eddy 
Simulations (LES) and Direct Numerical Simulations (DNS) for the computational 
analyses of high speed reacting flows. In the efforts related to LES, we have been 
primarily involved with assessing the performance of various modern methods based 
on the Probability Density Function (PDF) methods for providing closures for treating 
the subgrid fluctuation correlations of scalar quantities in reacting turbulent flows. In 
the work on DNS, we have been concentrating on understanding some of the relevant 
physics of compressible reacting flows by means of statistical analysis of the data gen- 
erated by DNS of such flows. In the research conducted in the second year of this 
program, our efforts have been focused on the modeling of homogeneous compressible 
turbulent flows by PDF methods, and on DNS of non-equilibrium reacting high speed 
mixing layers. Some preliminary work is also in progress on PDF modeling of shear 
flows, and also on LES of such flows. 

This report provides a summary of our achievements at the closing of the second 
year of this program. This research has been supported by NASA Langley Research 
Center under Grant NAG-1-1122. Dr. J. Philip Drummond, Theoretical Flow Physics 
Branch (TFPB), Mail Stop 156, Tel: 804-864-2298 is the Technical Monitor of this 
Grant 
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1 Introduction 


The need for the use of advanced computational methods for the analysis of high speed 
reacting flows is obvious [1, 2, 3]. Within the past thirty years, NASA Langley Research 
Center has been at the forefront in making use of advanced computational methods for the 
investigations of high speed turbulence, with and without the complex effects of chemical 
reactions. An overview of the work conducted at the Theoretical Flow Physics Branch reveals 
the extent of progress made by the scientists at NASA LaRC in this area of research. 

One of the recent contributions by NASA in the field of high speed combustion is due 
to Drummond [4, 5, 6]. This work, consisting of the first systematic study of high speed 
reacting shear flows, paved the way for the subsequent contributions made in this research 
area within the past few years. Since this early work, there have been extensive contributions 
made by other investigators on various aspects of high speed combustion. The extent of these 
achievements and the need for even further extensions have been so great that recently NASA 
in collaboration with ICASE, organized a workshop focused on identifying various means by 
which the area of high speed combustion can benefit from the current state of knowledge in 
the “simpler” field of low speed combustion. A proceeding of this workshop is available [7]. 
As asserted at the conclusion of this workshop, one of the major issues of crucial interest to 
NASA-LARC is the need for further developments and more substantiated work in utilizing 
advanced computational methods for the analysis of compressible turbulent combustion. 

A review of immediate needs in computational treatment of high speed flows is available 
[1]. According to this review, there is a need for further utilizations and improvements of 
the methods known as Direct Numerical Simulations (DNS) and Large Eddy Simulations 
(LES) for the computational treatment of compressible reacting turbulence. These two types 
of simulations have been labeled by this PI as “Model Free Simulations” [8], to emphasize 
the distinction between these methods and those based exclusively on “turbulence models” 
[9, 10]. Model free simulations have been the subject of wide utilizations by various groups 
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at NASA LaRC. This consists of contributions by both DNS and LES. Since the work of 
Drummond [4, 5], and the development of the computer code SPARK , the extent of achieve- 
ment has been promising. Some of the most recent contributions in this area are those in 
[11, 12, 13]. The SPARK code has been the subject of validations and improvements by 
various investigators at TFPB, NASA LaRC. Currently with the developments of modern 
numerical algorithms [14], advanced turbulence closures [15, 16], modern grid generations 
procedures [17], this code is being utilized as an efficient tool with accurate predictive capa- 
bilities for a variety of practical flow systems. This is clearly evident by the enthusiasm of 
the Aerospace industry in making use of the capabilities facilitated by this code [18]. 

There has also been a substantial interest at TFPB, NASA LaRC in developing and imple- 
menting LES for the treatment of compressible turbulent flows. This goers back to the early 
work of Erlebacher et al. [19], and since then LES have been the subject of wide utilizations 
by many scientist at this branch and also at ICASE [20, 21, 22]. Needless to indicate that the 
“Transition Group” at this branch is credited as the leader within the international scientific 
community in LES of compressible and transitional flows [23]. 

The outcome of earlier investigations and also those in progress clearly suggest that fore- 
seeable developments of advanced computational facilities will not be sufficient to relax 
the restriction of DNS to flows having small to moderate variations of the characteristic 
length and time scales. This is not an unknown problem, and has been well-documented in 
turbulence literature (for recent discussions see [24, 25, 8, 26]). The notorious limitations as- 
sociated with the direct simulation of high Reynolds number flows clearly indicates that the 
boundaries of applicability of DNS are, and will continue to be, somewhat restricted. Nev- 
ertheless, within its domain, DNS can be used to enhance our understanding of the physics 
of turbulent flows by (1) providing specific information concerning the detailed structures of 
the flow, and (2) providing a quantitative basis for evaluating the performance of turbulence 
closures. 

LES appear to provide a good alternative to DNS for computing flows having ranges of 
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parameters similar to those encountered in practical systems [27, 28, 29, 30, 31, 32, 33, 22, 
25, 8, 26]. The approach based on LES has a particular advantage over the turbulence 
modeling procedures in that only the effects of small-scale turbulence motion have to be 
modelled. Therefore, the construction of accurate “subgrid scale” closures is an important 
task in its implementation. The extensive experience gained during the past two decades in 
constructing turbulence models for the Reynolds- averaged equations of turbulent combustion 
have proven useful in the development of these closures. A major advantage offered by LES 
is that, subgrid closure modeling can be substantially simplified by performing computations 
over grids of different size [19, 21]. In this setting, the performance of a model on coarse 
grids can be directly evaluated by comparing its predictions with those obtained on fine 
grids. This procedure has been followed in many previous works [28, 19, 21, 20, 34] and the 
simulations, in many cases, have produced satisfactory subgrid closures (See [23] for a review 
of recent achievements). Most of these efforts have been towards constructing closures for 
non-reacting flows [30, 31], and it has been only recently that, the LES has been introduced 
for the analysis of reactive phenomena [34]. 

In the work conducted under this grant, we have continued our investigations on both LES 
and DNS related to compressible reacting flows. In the next section, the factors motivating 
research in the format proposed here are given, and in Section 3 a summary is provided of 
our accomplishments with this support at the conclusion of the second year of the research. 
The publications resulting from our efforts in Year 2 are listed in Section 3, together with 
the Awards and Honors received by the members of the research team supported by this 
Grant. 


2 Motivation 

The physics of high speed reactive flows is rich with many complexities [35], some of which 
the subject of intense investigations since the Apollo Project. A few examples of the physical 
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issues of current interest are the questions associated with the chemical and thermodynamics 
non-equilibrium effects, the cause and effect of turbulence, the interactions of turbulence 
and chemistry, the real gas effects at high temperatures, etc. This complexity necessitates 
investigations in many diverse areas consisting of all the branches of engineering. Having 
the discussions limited to technical issues, one of the most important elements of current 
interest is associated with the phenomenon of high speed turbulent combustion in internal 
flows [1]. Even in this relatively limited range, the physical complexities associated with 
“turbulence” is so broad that requires collaborative interactions amongst scientists within 
different disciplines in dealing with its cause and effects. 

An important tool in the analysis of high speed flows, as recognized by today’s Aerospace 
research community, is Computational Fluid Dynamics (CFD) [36, 37] (also see [38]). The 
power and the extent of capabilities of modern supercomputers have been instrumental in 
our efforts towards addressing and solving some of the existing problems in this field. 

As a result of collaborative efforts between SUNY-Buffalo and NASA-LaRC, we have been 
able to identify, and work on some important aspects of CFD utilizations in high speed 
combustion research. The particular aspect of the problem under investigation, in accord 
with our expertise, involves model free simulations of high speed turbulent combustion. 
Having the general guidelines and directions set by NASA-LaRC scientists, we initiated 
this work aimed at understanding several important issues in this direction. Our specific 
objectives from the start of this research has been to: (1) develop reliable mathematical 
closures for LES of high speed reacting flows, and (2) take advantage of the capabilities of 
DNS for further understanding of the physical issue associated with high speed transport. 
Our primary focus has been on the first task, while necessarily some efforts had to be devoted 
to the second task. 

In the next section, the summary of our to-date achievements are described. A more detailed 
descriptions of our achievements can be found in Appendix I through Appendix III. 
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3 Summary of Our Accomplishments to Date 


The subject matters discussed in this section are focused on both of the basic elements of 
the proposed research. Therefore, here they are discussed separately in order. 


3.1 LES and DNS of Homogeneous Reacting Turbulence 

It is now widely suggested that in developing a reliable subgrid closure for LES of turbulence, 
it is useful to consider a flow at the simplest possible condition. Namely, a homogeneous 
flow, with or without a spatial isotropy. This type of flow has been considered in almost 
all the previous works towards the developments of subgrid closures for LES [39, 40, 41, 
42, 19, 21, 20, 43]. The factors justifying this choice are based on the simple fact that the 
flow within the subgrid can be assumed isotropic regardless of the nature of the large scale 
behavior. Therefore, with the development of a turbulence model in a homogeneous setting, 
the generated results can be directly used in flows with strong spatial non-homogeneities. 
This line of reasoning is justified if the filtered portion of the transport variable constitute a 
sufficiently large portion of the “whole” flow, so that the assumption of spatial homogeneity 
within the un-filtered part is valid. 

In almost all previous work in developing subgrid closure for non-reacting flows, the approach 
based on “Moment Methods” has been adopted. The essential element of such methods is 
based on the revisions of the well-known Smagorinnsky [44] model (also see [45, 46]) which 
is somewhat analogous to the widely used gradient diffusion type closures in turbulence 
modeling [9]. The major modifications include the possible means by which the physical 
effects, such as compressibility, etc. are incorporated into the model. The extent of progress 
made to-date indicates that the state of knowledge necessary to incorporate such effects is 
not yet at a maturity level to justify absolute validity. For example, as indicated in [20], 
currently, only flows with low to moderate ranges of compressibility can be considered, and 
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the methodology is not capable of simulating flows at very high compressibility levels. 

Unfortunately, moment methods have not proven very effective in the modeling of turbulent 
reactive phenomena [47]. The results of extensive work within 1970’s-1980’s clearly reveals 
the shortcomings of such models primarily because of the lack of acceptable closures for all 
the ranges of the Damkohler number. Depending on the nature of the closure, a model which 
is reasonable in a “frozen flow” cannot be used for modeling of a chemically reacting flow 
under “equilibrium conditions,” and vise versa. This problem stems due to the drawbacks 
of the model to account for the effects of higher order unclosed moments at a selected level 
of truncation. 

The problem discussed above can be circumvented, to a large extent, by following the ap- 
proach based on the Probability Density Function (PDF) methods [48, 49, 50, 51, 8, 52]. 
The advantage of PDF methods in the modeling of reactive flows is due to the simple fact 
that with such models, all the contributions of the higher order moments are recovered by a 
single PDF transport equation. Therefore, once a model is developed there are no needs for 
constructing additional closures. The pleasing feature is that, this model can be developed 
in the simpler case of a non-reactive flow. Since the effects of chemical reactions appear in 
a closed form, the generated model can be subsequently used for the modeling of reactive 
phenomena [49, 51, 52, 53, 50, 8]. 

There are two general ways by which the PDF methods can be employed for solving turbulent 
combustion problems [10, 54, 8]: (1) Assuming the form of the PDF [48], (2) Solving a 
transport equation for the PDF [55, 56, 57]. Both of these approaches are well-known 
in combustion literature [15] and will not be elaborated upon here in detail. In the first 
approach, the PDF is parameterized with the knowledge of its first few (usually the first 
two) moments. This approach is obviously easier, but the moments must be provided by 
turbulence models. In the second approach, the consideration of the PDF transport equation 
is required. The problem with this approach, in addition to being more computationally 
demanding, is that some modeling is required to account for the the effects of the PDF 
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transport in the domain of the scalar composition. 

Knowing the general state of progress in PDF methods in turbulence modeling, we decided 
to examine both these PDF approaches in LES of turbulent reacting flows. That is, we 
initiated a parallel research program in which the potential capabilities and the limitations 
of both these methods were examined. In doing so, in many cases, we had to first focus 
on some of the properties of these models before being able to use them for the purpose of 
developing subgrid closures. Here, a summary is provided of all our findings in this aspect 
of our work: 

1. In the context of single-point PDF formulation, currently the Amplitude Mapping Closure 
(AMC) of Kraichnan [58, 59, 60, 61, 62, 63, 64] provides the most satisfactory results in the 
approach based on the PDF transport equations. This approach has proven more satisfactory 
than many of the previous models based on the so-called Coalescence/Dispersion (C/D) 
closure [65, 66, 50, 57, 67] in providing an acceptable evolution for the PDF transport. 

2. The main deficiency of the AMC (or any other single-point PDF model) is associated with 
its incapability in predicting the frequency (or length) scale of PDF evolution. The EDQNM 
(Eddy-Damped Quasi Normal Markovian) spectral closure [68] provides a reasonable means 
of overcoming this deficiency. In short, EDQNM provides a transport equation for the 
“covariance” of the scalar variable p( r) as a function of the “length scale.” The advantage of 
this closure is its capability of including the length scale information (through the parameter 
r) into the formulation. From this viewpoint, the closure is much more powerful tan the 
typical K — e [9] type closures. However, presently the use of the closure is limited to 
homogeneous and incompressible flows. 

3. The AMC is very effective in modeling of the mixing problem in homogeneous flows. 
Therefore, it is strongly recommended for statistical treatment of homogeneous reacting 
flows under chemically frozen and chemical equilibrium conditions. The use of the model is 
very convenient for the treatment of initially segregated reacting systems. In stoichiometric 
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mixtures, the limiting rate of fuel consumption in a chemical reaction of the type Fuel + 
Air — ► Products, the solution generated by AMC is of a simple algebraic form, and is very 
convenient for practical applications [48, 69, 70, 71]. However, it is not very useful for LES, 
because the reacting flow within the subgrid can hardly ever be assumed stoichiometric, 
albeit homogeneous [34]. The extension of the model to account for non-unity equivalence 
ratios has been made, but the final solution can be only expressed in terms of the integrals of 
“Degenerate Hypergeometric Functions.” In both cases, the predicted results compare very 
well with DNS results and are better than any other closure developed in the literature within 
the past forty years. Please see Appendix I for a detailed discussion of our achievements in 
this aspect of our work . 

4. The approach based on AMC requires a transport equation for the PDF. Therefore, its 
applications for LES requires numerical solution of an additional transport equation. It is 
shown that, in some circumstances, the family of probability frequencies of Pearson Type 
[72] provides a simple alternative to AMC. Namely, the /3-density of the first kind performs 
reasonably well in approximating the single point PDF of a conserved scalar property, if the 
evolution of its variance is known a priori. With the implementation of this density, it is 
also possible to develop a closed form algebraic solution for the rate of reactant conversion as 
indicated above. For a stoichiometric mixture the results are in terms of Gamma functions. 
For non-stoichiometric mixtures, the results can only be expressed in the form of Incomplete 
Beta Functions [73]. A mathematical justification has been provided of the similarity of the 
Pearson PDF’s with AMC generated frequencies. This justification fails if the PDF is not of 
a bimodal type. To the best of our knowledge, the extension of Pearson type distributions 
for the analysis of more complex distributions, is the subject of current research within the 
communities of statisticians and biometricians (see Appendix I). 

5. The use of the f3 density model has been extended for the modeling of non-homogeneous 
reacting flows. The simulated results are very encouraging provided that the evolution of 
the first two moments of the PDF are furnished externally. Again, the same limitations 
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indicated above, holds. Appendix II provides a complete description of our achievements in 
this aspect of our work. 

6. The Pearson family have been used for a priori analysis of subgrid modeling of the 
reacting flows in the limit of an infinitely fast chemical reaction. Our results suggest that 
in the absence of a better alternative, the beta density performs reasonably well in depicting 
the small scale behavior. However, as before, this does not advocate the generality of this 
density for practical applications. The results form this part of our work has been reported 
in [34], and has been included in our previous annual report to NASA LaRc (May 1991). 


3.2 DNS of High Speed Reacting Mixing Layers 

In the wake of our earlier works [12, 74, 13], we have continued our investigation of several 
issues of pertinence to the analysis of compressible reacting turbulence. This study has 
been motivated, to a large extent, in providing a computational complement to laboratory 
investigations of turbulent diffusion flames [75, 76]. We feel that, this is very important 
in view of the significant progress made within the past decade in the development and 
implementation of DNS for the analysis of unpremixed reacting flows [77, 78, 79, 80, 81, 82, 
83, 84, 85, 12] (see [86, 87, 1, 88, 89, 3] for recent reviews). The extent of progress made 
by these contributions has been very encouraging, thus justifying further utilization of the 
direct methods in the analysis of turbulent flames. 

As indicated before, in this aspect of our work, we were able to make use of present compu- 
tational capabilities provided by the SPARK code. In this period, we also became involved 
in utilizing other computer codes for the sole reason of making comparative assessments of 
the SPARK code, as well as enhancing our knowledge of the state-of-the-art in this sub- 
ject. Based on our work in the second year, we have been able to make the following solid 
conclusions: 
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1. With an idealized chemistry model of the type Fuel + Air — > Product + Heat , many 
fundamental issues in regard to the compositional structure of unpremixed spatially devel- 
oping reacting mixing layer can be investigated. More specifically, the results of our simula- 
tions do portray the same qualitative behavior as those of recent laboratory measurements 
[90, 91, 92, 93]. Appendix III provides a complete description of our achievements in this 
aspect of our work. 

2. For the first time, we have shown that exothermicity may actually act as to increase the 
rate of product formation in compressible reacting flows. This observation is not consistent 
with any of the previous DNS [94, 12] or experimental [95] results. The reason for this 
discrepancy is due to the effects of the Arrhenius kinetics model in our forced simulations 
and cannot be observed in the constant rate kinetics simulations and laboratory experiments 
[94, 12, 95]. This new observation, therefore, suggests that in the setting of a “turbulent” 
flame in which the “transitional effects” are not dominant, and where the chemical reaction 
is describable by an Arrhenius model, the heat release results in a higher reactant conversion, 
even though the rate of mixing may be somewhat less. Based on our observation, we recom- 
mend that the effects of exothermicity be assessed by means of laboratory measurements. 
These measurements must involve a reacting system whereby the rate of reaction conver- 
sion is temperature dependent (unlike that employed in [95]), and in which the large scale 
mixing intensity is not significantly affected by the heat release (see [94, 83, 88, 89, 84, 3] 
for discussions on the effects of heat release in low speed combustion, and [12, 74, 13] for 
discussions of such effects in high speed combustion). It is important to note that the issue 
would not be resolved by a mere comparison of the mixing characteristics between a reacting 
and a non-reacting layer. An appropriate analysis requires the consideration of two react- 
ing systems with the same hydrodynamic characteristics but with different magnitudes of 
the enthalpy of combustion. The extent of validity of our conclusion under more complex 
chemistry models can be determined by these experiments (see Appendix III). 

3. The result of our preliminary work in progress indicates that the concept of Laminar 
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Diffusion Flamelets [96, 97] does reasonably well in the modeling of low speed unpremixed 
reacting flows. That is the results generated by DNS of a turbulent non-equilibrium flame, 
scales well with the analytic solution of an opposed laminar jet configuration in relating the 
generated products to the instantaneous dissipations of the scalar field. However, it is con- 
cluded that in moderate chemical reactions, the “scatter” of the DNS results is substantially 
more than that to be completely described by the flamelet model. Please see Appendix III. 

4 Publications and Honors 

Within the second year of this research, two journal publications and one conference papers 
have resulted form this work. The publications are identified in the Reference list by numbers 
[98, 73]. The conference paper is identified by number [99]. 

Within the past year, one of the graduate students involved in this research, Mr. Steven 
Frankel, was named as the second place winner in the “Technical Paper Presentation” of the 
AIAA Northeastern Regional Conference. We would like to mention that the first winner in 
this part of the completion was Mr. Richard S. Miller, one of the other graduate students 
currently in our research group at SUNY-Buffalo. Also, in the “Technical Presentation” (no 
paper) part of this conference, Ms. Francine Battaglia (who also work with P. Givi), was 
named as the first prize winner. 
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5 Appendix I 


In this Appendix, we discuss the properties of the Amplitude Mapping Closure the Pearson 
Family of distributions for predicting the limiting rate of reactant conversion in homogeneous 
turbulent flows. The results provided in this appendix are useful in several regards: 

1. For the first time in forty year, closed form algebraic relations are obtained for predict- 
ing the maximum rate of reactant conversion in homogeneous turbulent flows. The major 
attraction of these relations is due to their applicability in both stoichiometric and non- 
stoichiometric (fuel-lean and/or fuel-rich) mixtures. As shown in this appendix, the results 
simplify substantially for unity equivalence ratio. The superiority of our results can be ap- 
preciated by a review of the most recent attempts in the modeling of homogeneous reacting 
turbulent flows (In Ref. [101], a review is presented of all the other closures available before 
this work), purpose). 

2. A direct consequence of our mathematical results is their applications for subgrid scale 
modeling in equilibrium reacting flows. As we have indicated in our previous Annual Reports 
to NASA LaRC, the flow within the subgrid is hardly ever stoichiometric. That is, the 
equivalence ratio within the subgrid is strongly dependent on the large scale mixing behavior. 
Therefore, it is very pleasing to have a closed form relation to predict the limiting rate of 
reactant conversion at all equivalence ratios. 

3. There is a careful discussion in this Appendix in extending this work for the modeling of 
non-equilibrium reacting flows with the implementation of multivariate statistical analyses. 

The extension of this work for predicting non- homogeneous flows are presented in Appendix 
II. 

The major part of this Appendix was prepared by Dr. Cyrus K. Madnia. 
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Abstract 


Closed form analytical expressions are obtained for predicting the limiting rate 
of reactant conversion in a binary reaction of the type F + rO — ► (1 + r) Product 
in unpremixed homogeneous turbulence. These relations are obtained by means of 
a single point Probability Density Function (PDF) method based on the Amplitude 
Mapping Closure (Kraichnan, 1989; Chen et a/., 1989; Pope, 1991). It is demonstrated 
that with this model, the maximum rate of the reactants’ decay can be conveniently 
expressed in terms of definite integrals of the Parabolic Cylinder Functions. For the 
cases with complete initial segregation, it is shown that the results agree very closely 
with those predicted by employing a Beta density of the first kind for an appropriately 
defined Shvab-Zeldovich scalar variable. With this assumption, the final results can 
also be expressed in terms of closed form analytical expressions which are based on 
the Incomplete Beta Functions. With both models, the dependence of the results on 
the stoichiometric coefficient and the equivalence ratio can be expressed in an explicit 
manner. For a stoichiometric mixture, the analytical results simplify significantly. 
In the mapping closure, these results are expressed in terms of simple trigonometric 
functions. For the Beta density model, they are in the form of Gamma Functions. In 
all the cases considered, the results are shown to agree well with data generated by 
Direct Numerical Simulations (DNS). Due to the simplicity of these expressions and 
because of nice mathematical features of the Parabolic Cylinder and the Incomplete 
Beta Functions, these models are recommended for estimating the limiting rate of 
reactant conversion in homogeneous reacting flows. These results also provide useful 
insights in assessing the extent of validity of turbulence closures in the modeling of 
unpremixed reacting flows. Some discussions are provided on the extension of the model 
for treating more complicated reacting systems including realistic kinetics schemes and 
multi-scalar mixing with finite rate chemical reactions in more complex configurations. 
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1 Introduction 


For the past forty years, since the early work of Hawthorne et al. (1949), estimation of 
the reactant conversion rate has been the subject of wide investigations in mathematical 
modeling of turbulent reacting flows. In unpremixed reacting systems, including diffusion 
flames, there are two factors by which this rate is influenced: (1) the speed at which the 
reactants are brought into the reaction zone, and (2) the rate at which they are converted to 
products through chemical reactions. These two mechanism are coupled through the mutual 
interactions of the fluid dynamics and the chemistry. The relative importance of the two 
mechanisms is characterized by the ratio of their appropriate timescales. This ratio is known 
as the Damkohler number ( Da ), and quantifies the characteristic frequency of the chemical 
reaction to that of the hydrodynamics. 

The role of the Damkohler number in the characterization of reacting flows is very important 
(Williams, 1985). As the magnitude of this number increases, the influence of chemical 
reactions become more pronounced and the reactants decay at a faster rate. This rate is 
bounded at the upper limit as the Damkohler number becomes infinitely large ( Da — ► oo). 
In this limit, the rate of reactant consumption is governed by the hydrodynamics, i.e. the 
reaction is “mixing controlled” and is determined by the speed at which the reactants axe 
brought into an infinitely thin reaction zone (Bilger, 1980). 

Obviously, with the assumption of an infinitely fast chemistry, it is not possible to account for 
many interesting issues associated with non-equilibrium effects in unpremixed flames (Libby 
and Williams, 1980). However, as indicated in the original pioneering work of Toor (1962), 
and later by O’Brien (1971); Bilger (1980), it is very important to have a prior estimate of 
this limiting rate in practical modeling of reacting flow systems. This is simply due to the 
fact that in this limit, the problem reduces to the simpler problem of “mixing” in which its 
analysis is far less complicated than that of an equivalent reacting system. Furthermore, it 
provides an accurate estimate for those reacting systems in which the chemical reaction time 
is much smaller than the characteristic time associated with the hydrodynamics. 
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In most analysis of turbulent reacting flows, the “statistical mean” value of the reactant 
conversion rate is of practical importance, and development of an appropriate turbulence 
model which can predict this mean rate has been the subject of extensive investigations 
(for reviews see Toor (1975); Brodkey (1981); Libby and Williams (1980); Williams (1985)). 
Amongst the theoretical tools developed, it is now firmly accepted that the approach based on 
the single-point Probability Density Function (PDF) of the scalar quantities is particularly 
useful, and this approach has been very popular in modeling of the reactant conversion 
in a variety of turbulent reacting flow systems (Pope, 1991; Kollmann, 1990; Pope, 1985; 
Pope, 1990). The advantage of PDF methods is due to its inherent capability to include 
information for all the statistical correlations amongst the scalar field. Therefore, once the 
PDF (or the joint PDF) of the scalar variables is determined, all the relevant statistics of 
the field are available without a need for additional closures. 

The most logical and systematic means of determining the PDF involves the solution of 
an appropriate transport equation governing its evolution in a given reacting flow. In this 
equation, due to the nature of the formulation, the effects of chemical reaction appear in 
a closed form (Pope, 1976), regardless of its degree of complexity. However, the influences 
of molecular action cannot be fully described, and can be treated only by means of em- 
ploying an appropriate closure. As noted by Pope (1991), in most previous applications, 
this problem has been circumvented through the use of the Coalescence/Dispersion (C/D) 
models. Examples of such models are the early C/D prototype of Curl (1963), the Lin- 
ear Mean Square Estimation (LMSE) theory (O’Brien, 1980), the closure of Janicka et ai 
(1979), among others (Pope, 1982; Pope, 1985; Kosaly and Givi, 1987; Givi, 1989). De- 
spite their advantageous characteristics, the shortcomings associated with the C/D closures 
in the probabilistic description of scalar transport are well recognized. Namely, none of 
the aforementioned models predict an asymptotic Gaussian distribution for the PDF of a 
conserved scalar variable in homogeneous turbulence (Pope, 1982; Kosaly and Givi, 1987; 
Pope, 1991); and those which are capable of doing so ( e.g . Pope (1982)), do not predict the 
initial stages of mixing correctly (Kosaly, 1986). 

Recent development of the Amplitude Mapping Closure by Kraichnan and co-workers (Kraich- 
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nan, 1989; Chen et a/., 1989) (see also Pope (1991)) has provided a promising way of alleviat- 
ing some of the problems associated with the C/D closures. This closure, in essence, provides 
a means of accounting for the transport of the PDF in composition space, and its validity 
and physical applicability have been evidenced in a number of comparisons against data gen- 
erated by means of both direct numerical simulations (Pope, 1990; Pope, 1991; Gao, 1991a; 
Gao, 1991b; Madnia et al., 1991b; Jiang et al., 1992), and laboratory experiments (Frankel 
et al., 1992). These results suggest that, at least in the setting of an isotropic turbulent flow, 
this closure has some superior features over all the previous C/D type models. 

Based on this demonstrated superiority, our objective here is to further examine the proper- 
ties of this closure and to assess its capabilities for applications in modeling of unpremixed 
turbulent reacting flow systems. In particular, it is intended to provide a reasonably simple 
recipe that can be used in conjunction with this closure for predicting the limiting rate of 
mean reactant conversion. However, since this is the first study of this type, and due to 
mathematical complexities (that soon become apparent), we have made some simplifying 
assumptions, which are indicated here at the outset. Firstly, we consider an idealized ir- 
reversible binary reaction of the type F + rO — * (1 + r) Product with initially segregated 
reactants ( F and O ). In accordance with the discussions above, only the maximum rate 
of reactant conversion is considered, implying that the magnitude of the Damkoher num- 
ber is infinitely large. Secondly, the turbulence field is assumed statistically homogeneous. 
Thirdly, all the chemical species are assumed to have identical and constant thermodynamic 
properties. Finally, the flow field is assumed isothermal in which the dynamic role of the 
chemical reaction on the hydrodynamic field is ignored. 

With all these assumptions, the reacting system considered is obviously an idealized proto- 
type of conventional combustion systems. However, it does provide a good model for dilute 
reacting systems in typical mixing controlled plug flow reactors (Toor, 1962; Bilger, 1980; 
Toor, 1975; Hill, 1976; Brodkey, 1981). Moreover, because of the mathematical complexi- 
ties, even in this simple case, it is deemed necessary to analyze this simplified system before 
considering more complex scenarios. Nevertheless, the model is capable of accounting for 
arbitrary values of the stoichiometric coefficients and for any equivalence ratio. This allows 
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capturing of many interesting features, as will be demonstrated. Some issues regarding fu- 
ture extension of this work will be discussed. At this point, it is sufficient to indicate that the 
approach followed here can be directly extended for treating more complicated flows with 
relaxation of all of the assumptions made here. However, in those cases, the final results 
cannot be generally expressed in terms of simple mathematical relations, and require numer- 
ical solution of the governing equations. Pope (1991) describes the procedure by which the 
transport equations generated by the mapping closure can be treated numerically. 

For the idealized case of initially segregated reactants F and 0, as indicated above, the 
initial marginal PDF’s of their mass fractions are composed of “delta functions”. Therefore, 
it is also speculated that the approach based on an assumed probability distribution may 
also provide a reasonably good closure. Therefore, in addition to the mapping closure, a 
member of the family of Pearson type frequencies is also considered. The results obtained 
by this frequency axe compared with those of the mapping closure, and are also assessed 
against data generated by means of DNS. 

In the next section, the problem under consideration is outlined along with the mathematical 
basis by which the single- point PDF methods are used. In Section 2.1, the salient features 
of the mapping closure at the single-point level are analyzed including a discussion on the 
formalities of the closure for our purpose. With this closure, a closed form analytical ex- 
pression is provided for the limiting bound of the reactant conversion. This limiting rate is 
also predicted by means of the Beta density model in Section 2.2. In both these subsections, 
the simplifications for the case of a stoichiometric mixture are made, motivating the use of 
our final “simple” analytical expressions in practical modeling of stoichiometric plug flow 
reactors. In Section 3.1 the results predicted by both models are compared against those 
generated by DNS of a three-dimensional homogeneous turbulent flow. In Section 3.2 some 
discussions are presented highlighting the implications of these results in turbulence model- 
ing. This paper is drawn to a close in Section 4 with some discussions for possible future 
extension of the two models for the statistical description of more complicated chemically 
reacting turbulent flows. 



2 Formulation 


With the assumptions described in the introduction, the statistical behavior of the reacting 
field in the reaction F+rO —* (1 +r)Product can be related to the statistics of an appropriate 
conserved Shvab-Zeldovich mixture fraction, J (Bilger, 1980). This mixture fraction can be 
normalized in such a way to yield values of unity in the fuel F stream and zero in the oxidizer 
O stream. For the purpose of statistical treatment, we define V F ({, 'Po{£, t) and Vj(£, <), 

respectively, as the margined PDF’s of the mass fraction of F, the mass fraction of O, and the 
Shvab-Zeldovich variable J . For initially segregated reactants with no fuel in the oxidizer 
stream (and vice versa), the initial conditions for the margined PDF’s of the mass fractions 
of the two reactants are given by: 


Pf(Z, 0) = W F S{i - Fi) + W o 6(0, 


PoU, 0 ) = Wo6(t - Oi) + W f 6(0, o < * < 1. (1) 

Here, Fi and Oi denote the initial mass fractions of the two species in the two feeds, and W F 
and Wo represent the relative weights of the reactants at the initial time (i.e. the area ratios 
at the inlet of a plug flow reactor). With the normalized value of the mass fractions equal to 
unity at the feeds, i.e. Fi = 0, = 1, the stoichiometric value of the Shvab-Zeldovich variable, 
is determined from the parameter r. With the assumption of an infinitely fast chemistry, 
the marginal PDF’s of the two reactants are related to the frequency of the Shvab-Zeldovich 
variable (Bilger, 1980; Kosaly and Givi, 1987): 




Vo(U) = ( 2 ) 


6 


Here, £ > 0, and 


<v(o= 

J — OO 


M‘) = rvjd.tw = 1-<M*)- (3) 

JJ. I 

The initial condition for the PDF of the Shvab-Zeldovich variable is given by: 

iV«,o) = w^K-i) + W«). (4) 

Equation (4) implies that < J > (f = 0) = Wp. Since J is a conserved variable, its mean 
value remains constant, i.e. < J > (t) =< J > (0) = W F . The integration of Eq. (2) yields 
the temporal variation of the statistics of the species field at all times, if the PDF’s of J are 
known. As indicated above, in the setting of a mixing controlled reaction this PDF provides 
all the desired statistical properties of the reacting field. 


2.1 Amplitude Mapping Closure 

The implementation of the amplitude mapping closure involves a mapping of the random field 
of interest £ to a stationary Gaussian reference field <f> o, via a transformation £ = x(^o ? 0- 
Once this relation is established, the PDF of the random variable £, 'P(^), is related to 
that of a Gaussian distribution. In homogeneous turbulence, the transport equation for this 
function satisfies (Chen et al., 1989; Pope, 1991): 


dx , d 2 X 
dr~ dft 


( 5 ) 


In this equation, r is a normalized time within which the scalar length scale information 
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is imbedded. The relations between this time and the physical time, i.e. r(t), cannot be 
determined in the context of single-point PDF description and must be provided by external 
means (Pope, 1990; Jiang et al., 1992). For the case considered here, with the initial PDF 
of the variable J given by Eq. (3), the corresponding form of the initial mapping is: 


x(<£o,0) = H(<t>o - <ff), — oo < <f> 0 < oo, (6) 

where H is the Heaviside function, and 4> m is a measure of the initial asymmetry of the initial 
PDF around the ensemble mean of the variable J\ 


<f>* = y/2ed~ l (l -2 < J >), (7) 

where “erf” denotes the error function. The mapping function is obtained by solving Eq. 
(5) subject to initial condition (6). The general solution of this equation has the form (Gao, 
1991a): 


\[<P + 1 r°° 


[OO 

(<^oe T -y) 2 (l + Q 2 ) 

/ XU/> u J ex P 

r — oo 

2Q 2 


dy. 


( 8 ) 


where Q{t) = ^/exp(2r) - 1. Inserting Eq. (6) for x(^o,0) in Eq. (8), we have: 


x(<t>o, r) = ^ (1 + erf(a^ 0 + *>)] , 


( 9 ) 


where, 


a(r) 


1 u ^ -fv'i+g* 
v/2S’ 6(T) “ V2S 


( 10 ) 


Finally, the solution for the PDF of J is determined directly from the mapping relation 
between the physical field £ and the Gaussian reference field <f>o: 
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( 11 ) 


Fj (x(^o» r ). r ) = Pa(to) (Jj^j • 

Here, Pg denotes the PDF of a standardized Gaussian distribution, i.e., Pg{4> o) = exp(— <t>l/2). 

A combination of Eq. (11) and Eq. (9), yields the final result for PDF of the Shvab-Zeldovich 
variable: 


Fj (x(^o> r )> r ) = ^exp 


(<^ 0 g- r - r? 

2(1 — e _2r ) 



( 12 ) 


With a combination of Eqs. (12) and (2), all the pertinent single-point statistics of the 
reacting field are determined. The most important of these statistics are the ensemble 
mean values of the reactants’ concentrations. These mean values are obtained directly by 
integrating their respective PDF’s. The intermediate steps in deriving these relations are not 
presented but are provided by Frankel (1992). Here, only the essential steps are presented. 
For the mean fuel concentration, < F > the first part of Eq. (2) reads: 


<F>(t)= l 1 F({)Vj(Z,T)d( = r F(x,r)P G (<f>o)d4>o, (13) 

Jj. , ) 

where the lower limit of the last integral corresponds to the value of <f > 0 at which \ is equal to 
the stoichiometric value of the Shvab-Zeldovich variable. Evaluating this limit from Eq. (9), 
Eq. (13) can be analytically integrated. This is possible by representing the error function in 
the form of its definition, and performing the resulting definite, double exponential integral. 
The results after extensive algebraic manipulations yield: 


< F > (r) = 


(1 - 2 J„) 

4(1 -J„) 



+ tv'm! - j .) exp 



2 



( 14 ) 


Similar expressions can be obtained for the mean oxidizer concentration: 
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<0>(t) 2 (1 2J J 


. l - erf (^f f )]-^: exp (“^ _ ^ _ 9 e - (15) 


H 




In these equations, c is related to the stoichiometric coefficient, 


c = — erf- 1 (2J gt - I) 
a 


(16) 


and 



e ± = /Vx + | 



ac 

VZT2 


X>-i 



V2 = y +^, 


p 


y± = ^2 acy 2 


b 


c 

a 


(17) 


Here, X>_ 2 and D_ x are, respectively, the Parabolic Cylinder Functions of order —2 and -1, 
belonging to the family of degenerate hypergeometric functions (Abramowitz and Stegun, 
1972). 


Due to nice properties of the degenerate hypergeometric functions, many of the interesting 

features of Eqs. (14)-(15) can be depicted. First, simple manipulation of this equation shows 

that for a stoichiometric mixture, < J >= J st , both reactants decay at the same rate, i.e. 

<£>!*! — = and for a non-stoichiometric initial condition, the limits of the 

<F>(0)=VV> <O>(0)=W o ’ 

mass fractions as r, Q — ► oo, asymptote to: 


f° 

Lim^-cc < F>(t)= \ 

l i-J« 


Jst > <J> 

Jst < <J > 
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(18) 


*r-» 


Linv.c-oo < O > (r) = 



< < J > 

J a t > <J > 


These limits are obtained by employing the Taylor series expansion of the relevant functions 
as Q — ► oo, and indicate the limiting bound of the unconsumed reactants in both fuel rich 
and fuel lean mixtures. While these limiting behaviors are rather trivial from a physical 
standpoint, in a computational procedure one must make sure that they are satisfied. Be- 
cause of the mathematical properties of the Parabolic Cylinder Functions, these limiting 
cases can be realized in our computational procedure in a relatively easy manner. It would 
be very difficult to obtain these limiting behavior numerically in an integration procedure 
within the original unbounded domain. 

At .first glance, Eqs. ( 14)-( 15) may appear somewhat complicated. However, due to nice 
mathematical properties of the Parabolic Cylinder Functions (Abramowitz and Stegun, 
1972), these equations can be integrated rather easily within the finite domain (0 < y < 1). 
For fuel-lean or fuel-rich mixtures, the integration can be only done by means of employing 
numerical methods. However, for a stoichiometric mixture, the results simplify further as 
demonstrated below. 



Stoichiometric Mixture: 

For practical applications in stoichiometric plug flow reactors, the equations simplify con- 
siderably. For a stoichiometric mixture, and an initially symmetric PDF around the mean 
value, (i.e. < J >= = 1/2), both parameters b and c are zero. Under this condition, the 

first terms on RHS of Eqs. (14)-(15) drop. Knowing £L 2 ( 0) = X>-i(0) = 1, the remaining 
terms yield: 


< F > (t) _ < 0 > (t) _ 1 f l dy _ ^ _ 2arctang(r) 

< F > (0) _ < 0 > (0) " Tsfia Jo 2y 2 + ^ i r 


( 19 ) 
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The simplicity of this equation is noteworthy and very pleasing. Because of this simplicity, 
Eq. (19) is strongly recommended for engineering predictions of the reactant conversion rate 
in stoichiometric homogeneous flows, such as the plug flow reactors considered in numerous 
previous investigations (Toor, 1962; Toor, 1975; Brodkey, 1981; Hill, 1976; 0 Brien, 1971; 
Kosaly and Givi, 1987). 


2.2 Beta Density Model 

For initially segregated reactants, the initial PDF of the Shvab-Zeldovich variable is com- 
posed of two delta functions at the extreme limits of the variable. Therefore, it is proposed 
that the family of Pearson (1895) frequencies may provide a reasonable means of estimating 
this distribution at all times (Frankel et a/., 1991; Girimaji, 1991a). The appropriate form 
of the Pearson distribution is in this case the Beta density of the first kind. This density has 
been employed in the statistical description of turbulent reacting flows by (Rhodes, 1975; 
Jones and Priddin, 1978; Lockwood and Moneib, 1980; Peters, 1984; Janicka and Peters, 
1982) amongst others (for recent reviews, see (Givi, 1989; Priddin, 1991)). For an ini- 
tially non-symmetric PDF, the Beta density corresponds to the Pearson Type I and for the 
symmetric case to Pearson Type II. The relevance of the latter in modeling of molecular 
mixing from an initial symmetric binary state has been described by Madnia et al. (1991a); 
Girimaji (1991a). In both cases, the PDF of the Shvab-Zeldovich variable is represented by 
(Abramowitz and Stegun, 1972): 


where B((3i,/3 2 ) denotes the Beta function, and the parameters & and /3 2 are dependent on 
the mean and the variance of the random variable 3 . In applications to the mixing controlled 
reaction considered here, we assume that the PDF of the Shvab-Zeldovich variable always 
retains a Beta distribution. Thus all the statistics of the reacting scalar are subsequently 
determined. The ensemble mean values are determined by a combination of Eqs. (20) and 
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(2). Following the same procedure as that described in the section on the mapping closure, 
after some manipulations the final results can be expressed as (Frankel, 1992): 


' > <» - (srs - 4 “ ' - 1 ™ 


(A + fli) 


< 0 > W = t£t + ( x - (A + A)J„) 2> “ (A ’ A) ’ 


( 22 ) 


where J denotes the Incomplete Beta Function (Abramowitz and Stegun, 1972). 

Due to nice mathematical properties of the Beta Function, the final results are cast in 
terms of its integral. In this case, however, the integral can be expressed in terms of the 
known special Incomplete Beta function. The mathematical properties of this function are 
well known, and the expressions above are conveniently amenable to numerical integration 
(Frankel, 1992). Again, the physical limiting conditions discussed before are realized by Eqs. 
(21)-(22). That is, in a stoichiometric mixture, both reactants decay at the same rate; and 
in lean or rich mixtures, the same limiting conditions as those in Eq. (18) are realized. 


Stoichiometric Mixture 

Again, in the case of an initially symmetric PDF under stoichiometric conditions, the final 
expressions become simpler. Under this condition, fli = (3 2 , and knowing T x / 2 (x,i) = 1/2, 
Eqs. (21)- (22) reduce to: 

<F>{t) _ < 0 > (t) 1 r(g) (23) 

< F > (0) < 0 > (0) v^ r (0 + D’ 

where g is half the inverse of the normalized variance of the Shvab-Zeldovich variable, and 
T denotes the Gamma function. 
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3 Results 


The final forms of Eqs. (14), (15), (19), (21), (22), and (23) are gratifying since they provide a 
relatively simple and effective means of estimating the maximum rate of reactant conversion 
in homogeneous reacting flows. As indicated before, the mathematical operations leading to 
these equations are somewhat involved, but the final results can be conveniently expressed 
in terms of known special functions. However, since both models are based on single-point 
PDF descriptions, these equations are not in a complete closed form and are dependent on 
the parameters g and/or Q. In this context, this parameter cannot be determined by the 
model and must therefore be specified by external means (Pope, 1991; Jiang et a/., 1992; 
Frankel et a/., 1992). This deficiency is not particular to the two models considered here, 
and exists in any single-point statistical description, including all of those based on the C/D 
models. 

The extent of validity of these simple relations can be demonstrated by a comparison between 
the model predictions and the data obtained by means of DNS. The comparison is made here 
for several values of < 3 > and J7jt for the purpose of demonstration. In this comparison, 
the magnitudes of the normalized variance of the Shvab-Zeldovich variable are matched with 
those of DNS. This implies that for given values of < J > and J7^ t , the parameters g , G,0i, 
and & are provided externally from the DNS data. With this provision, the model prediction 
results can be directly assessed against DNS data. 

The DNS procedure is similar to that of previous simulations of this type. For a detailed 
description, we refer the reader to Madnia and Givi (1992). The subject of the present DNS 
is a three-dimensional periodic homogeneous box flow under the influence of a binary reaction 
of the type described above. The initial species field is assumed to be composed of out of 
phase square waves for the two reactants F and O. The computational package is based on 
the modification of a spectral-collocation procedure using Fourier basis functions developed 
by Erlebacher et al. (1990a) (also see Erlebacher et al. (1987); Erlebacher et al. (1990b)). 
The hydrodynamic field is assumed isotropic, and is initialized in a similar manner to that 
of Erlebacher et al. (1990a); Passot and Pouquet (1987). The turbulent field is of decaying 
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nature in that there is no artificial forcing mechanism to feed energy to low wave numbers. 
The code is capable of simulating flows with different levels of compressibility (Hussaini et 
ai, 1990). Here, only the results obtained for a low compressible case are discussed, since 
most previous analyses of plug flow reactors have dealt primarily with incompressible flows 
(Toor, 1975; Hill, 1976; Brodkey, 1981; Leonard and Hill, 1988b; Leonard and Hill, 1991). 
The resolution consists of 96 collocation points in each direction. Therefore, at each time 
step 96 3 is the sample size for statistical analyses. With this resolution, simulations with a 
Reynolds number (based on the Taylor microscale) of Re\ « 41 are attainable. The value of 
the molecular Schmidt number is set equal to unity. 


3.1 Validations 

The statistical behavior of the scalar field is depicted by examining the evolution of the 
PDF’s of the Shvab-Zeldovich variable J. These are shown in Figs. 1 and 2 at times close 
to the initial (£1) and the final (<2) states. These figures correspond, respectively, to the 
cases of an initially symmetric (< J >= 1/2) and non-symmetric (< J >= 0.625) PDF’s. 
At the initial time, the PDF is approximately composed of two delta functions at J = 0, 1 
indicative of the two initially segregated reactants, F and O. At later times, the PDF evolves 
through an inverse-like diffusion in composition space. The heights of the delta functions 
decrease, and the PDF is redistributed at other J values in the range [0, 1], and subsequently 
becomes centralized around the mean value. Proceeding further in time results in a sharper 
peak at this mean value, and in both cases, the PDF can be approximated by a Gaussian 
distribution. The trend for the symmetric case is the same as that presented in earlier DNS 
studies (Eswaran and Pope, 1988; Givi and McMurtry, 1988). For the non-symmetric case, 
there are no DNS data in the literature, but the present results verify that the asymptotic 
PDF can still be approximated by a Gaussian distribution. 

The PDF’s obtained by the mapping closure and those assumed by a Beta density are 
also presented in Figs. 1-2. In these figures, the model PDF’s are parameterized with the 
same first two moments obtained from DNS. In this parameterization, only the normalized 


magnitude of the variance of the models are forced equal to that of the DNS and no attempt 
was made to account for the departure from the “exact” initial double delta distribution in 
DNS. With this matching, nevertheless, the results clearly indicate that the model predictions 
compare very well with the DNS results. Also, both models yield an asymptotic Gaussian 
PDF. 

The temporal variation of the ensemble mean of the reactants’ mass fractions by the two 
models are compared against those of DNS in Figs. 3-4. These figures correspond to the 
two cases of symmetric and non-symmetric initial PDF’s, respectively. In the symmetric 
case, under stoichiometric conditions, the results are simply obtained from the analytical 
expressions in Eqs. (19) and (23). In the non-symmetric case, the numerical integration of 
Eqs. (14)-(15), and the evaluation of the Incomplete Beta Function in Eqs. (21 )-(22) are 
necessary. In both cases, the agreement of the models with the DNS data is noteworthy. 
Also, a comparison between parts (a) and (b) of Fig. 4 shows that as the magnitude of t 
decreases, the rate of oxidizer consumption becomes faster. 

The agreement noted above is pleasing but not very surprising. It follows from the compati- 
bility of the model PDF’s and those of DNS, at least for the case of binary mixing considered 
here. This finding is not new and has been well documented in previous works, at least those 
considering an initially symmetric PDF (Pope, 1990; Pope, 1991; Madnia and Givi, 1992; 
Madnia et a/., 1991a; Girimaji, 1991a). However, a nice feature of the models is the explicit 
form of the final equations expressing these statistical quantities. Supported by this quan- 
titative agreement, it is proposed that in the absence of a better alternative, the relations 
obtained above to be used as an explicit simple means for predicting the maximum rate of 
reactant conversion in homogeneous reacting systems. 

Despite the simplicity of these equations and their ease of applications, it must be mentioned 
that these equations predict the rate of mean decay of the reactant far better than all of the 
previous turbulence closures based on the C/D models (see Givi (1989)). This is particularly 
advantageous in that this evaluation can be made by a simple algebraic relation, whereas 
the C/D implementations usually require more expensive numerical simulations. Even for 


16 


non-unity equivalence ratios, the numerical integration required by the two above models is 
considerably less computationally demanding than those of the C/D models. The only input 
in these models, similar to those in C/D closures, is the variance of the Shvab-Zeldovich 
variable. This is provided here by means of DNS. In an actual implementation, this variance 
can be obtained from experimented data or by means of an appropriate turbulence model 
(Frankel et al., 1992). The provision of such data is not very difficult since they can be 
obtained in the setting of a non-reacting flow. 


3.2 Applications 

The relations obtained here can be used in determining the extent of validity of other con- 
ventional closures for predicting the limiting rate of reactant conversion in turbulent flows. 
As an example, here we consider the model based on the famous hypothesis of Toor (1962); 
Toor (1975), which has received considerable attention in practical modeling of unpremixed 
homogeneous reacting systems (Bilger, 1980; Brodkey, 1981; Leonard and Hill, 1987; Leonard 
and Hill, 1988a; Leonard and Hill, 1988b; Kosaly and Givi, 1987; Kosaly, 1987; Givi and 
McMurtry, 1988; McMurtry and Givi, 1989; Givi, 1989) . According to this hypothesis, in an 
isothermal homogeneous reacting turbulent flow, the decay of the unmixedness denoted by 
< a'V >, where the prime quantities indicate fluctuation from the ensemble mean value, is 
independent of the magnitude of the Damkohler number. This implies that the normalized 
unmixedness parameter, defined by: 


_ < o'fe' > (0|p, 

< a'V > {t)\ Da = 0 ' K ] 

is the same under both reacting and non-reacting conditions, i.e. Z(t) = constant = 1 for 
all values of Da. In previous DNS assessments of this hypothesis, it has been shown that for 
the case of initially segregated reactants this model cannot be employed, and the normalized 
unmixedness ratio depends on the nature of mixing and the magnitude of the Damkoher 
number (Givi and McMurtry, 1988; McMurtry and Givi, 1989). In particular, it has been 
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demonstrated that even for Da — *• oo, while the normalized unmixedness is equal to unity 
at the initial time, its limiting lower bound depends on the asymptotic frequency of the 
Shvab-Zeldovich variable. That is, Z(t = 0) = 1 > Z(t) > Z(t — » oo) = C, where C is 
the lower limiting bound. For an asymptotic Gaussian distribution, it can be easily shown 
that for a mixture under stoichiometric conditions, the lower bound limit asymptotes to the 
constant value C = 2/ ir (Kosaly, 1987; Givi and McMurtry, 1988). 

This deviation from unity can be realized by the two models considered. With the mapping 
closure, under symmetric stoichiometric conditions, from Eqs. (2), (14), and (15), following 
the same integration procedure as before, it is shown that: 


< a!V > {Da -» oo) _ 2 (arctan(^)) 2 

< a'hf > (Da = 0) * arC tan ( g ^ +2 -) 


(25) 


For the Beta density model, the corresponding form of the unmixedness ratio for /3i = ft? is 
given by: 


Zff- 



(26) 


In Eqs. (25)-(26), the subscripts m, /? are added to denote the mapping and the Beta density 
model. It is easy to show that these equations satisfy the correct limiting conditions for a 
stoichiometric mixture. This is for both the initial time, i.e. the inlet of the reactor, 


Lim[$_o].Z m — Limj^^ij^ij 1, 


(27) 


and at large distances form it: 


Limp^oo] Z m 


Lim[j— -°°] 2p 



(28) 
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The latter limiting condition cannot be realized in any of the previously used C/D models, 
or by means of the Toor ’s models (McMurtry and Givi, 1989; Givi, 1989). 

The results based on the applications of Toor’s model become less accurate for non-stoichiometric 
mixtures. For equivalence ratios other than unity, with the depletion of one of the reactants, 
the unmixedness parameter approaches zero faster. This is demonstrated in Fig. 5 for several 
values of the equivalence ratio (7). Note that as the magnitude of this ratio increases above 
one, the unmixedness ratio goes to zero more rapidly. For an unity equivalence ratio, the 
correct asymptotic value of 2/x is realized. 


4 Extensions for Modeling of More Complex React- 
ing Turbulent Flows 

Despite the pleasing features of our simple mathematical expressions, there are several re- 
stricting assumptions which were necessarily imposed in deriving these equations. Here, we 
would like to address the ramifications associated with these assumptions, and to provide 
the means of overcoming them in future extensions of these two models. 

Firstly, due to the assumption of infinitely fast chemistry, only the maximum rate of the re- 
actant conversion is obtained. While this rate is of crucial interest in describing unpremixed 
flames, from both a theoretical standpoint and for practical applications (Givi, 1989; Mc- 
Murtry and Givi, 1989; Kosaly and Givi, 1987; Toor, 1962; Toor, 1975; O’Brien, 1971; 
Bilger, 1980; Williams, 1985), the model is not capable of describing some of the important 
features of the turbulent flames, especially those under non-equilibrium conditions. The ex- 
tensions to finite rate chemistry, reversible reactions, non-equilibrium flames and multi-step 
kinetics systems require numerical integration of the PDF transport equation. For these 
cases, the problem cannot be mathematically reduced to that of keeping track of a single 
scalar variable (like J), and requires the use of multivariate statistical descriptions. For this, 
the implementation of mapping closure is relatively straightforward since it provides a trans- 
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port equation for the joint PDF’s of the scalar variable in a multivariable sense (Pope, 1991; 
Gao and O’Brien, 1991). However, it is not presently clear how to devise an efficient compu- 
tational procedure, typically based on Monte Carlo methods (Pope, 1981), for the numerical 
treatment of these equations. Some work in this regard is currently under way (Valiho and 
Gao, 1991; Valino ei a/., 1991). 

The extension of assumed distributions based on the Beta density for treating multi-scalars 
is more straightforward but less trivial to justify. The most obvious means is to implement 
the multivariate form of the Pearson distributions. The direct analogue of the Beta density 
is the Dirichlet frequency (Johnson, 1987; Narumi, 1923) . For a mixture composed of N + 1 
species, the joint PDF of N mass fractions (rp i, 02> • ■ • *Pn) is described in terms of a N- variate 
density of the form: 


Pty 1,^2,- ■■ i W 


r(/?i + th + : 'j£n± i) 

r(A)r(^ 2 ),...r(^ + i) W2 


. . . (i -ih - * - . . . M 0 ”*' 

(29) 


subject to the physical constraint: 


JV+l 

E v>,- = i ( 3 °) 

i=i 

The application of this density in modeling of multi-species reactions has been nicely dis- 
cussed by Girimaji (1991a). Due to the mathematical properties of the Gamma function, 
this density is pleasing from a mathematical viewpoint and most statistical cross correla- 
tions of the random variables (ipi,ip 2 , • • ■) can be conveniently obtained by means of simple 
analytical relations (Frankel, 1992). Some points in this regard has been made by Girimaji 
(1991b). However, the use of the Dirichlet frequency cannot be justified for modeling of un- 
premixed reacting flow in a general sense. In fact, it can be shown that in the case of a binary 
irreversible reaction of the type considered here, if the joint PDF of the two species can be 
assumed to be of a bivariate Dirichlet then the univariate PDF of the Shvab-Zeldovich vari- 
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able cannot be approximated by a Beta density. Moreover, there is no way of implementing 
this density directly for modeling of non-equilibrium flames, involving strong temperature 
variations. This is simply due to the additivity constraints of this density requiring the unity 
sum of the normalized random variables (Eq. 30). 

Secondly, the mathematical derivations presented here are only valid for initially segregated 
reactants. In both models, the complete segregation facilitates significant simplifications 
of the final equations. This assumption is not unrealistic and is compatible with that in 
majority of previous works on unpremixed reacting flows (Toor, 1975; Brodkey, 1981; Bilger, 
1980). For other initial conditions, e.g. partial premixing of the reactants, or non-delta 
like distributions, one must resort to numerical integration of the PDF transport equation. 
Again, an appropriately devised numerical procedure can accommodate such conditions. But 
the use of a Beta density (or any other assumed distributions) cannot be justified for other 
complex initial conditions. 

Thirdly, the final mathematical expressions are only valid in the setting of a homogeneous 
flow. The extension to inhomogeneous flow predictions is also straightforward, but requires 
numerical integration of the modelled equations. Both models can be directly implemented 
into appropriately devised numerical procedures. The mapping closure can be invoked in 
the mixing step of a fractional stepping procedure, similar to that of typical Monte Carlo 
procedures (Pope, 1981). The Beta density requires modelled transport equations for the 
low order statistics of the reacting field. These equations include the required information 
pertaining to the spatial inhomogeneity of the flow, through the parameters /?i , ^ 2 » • • With 
these informations, all the higher order statistics of the reacting field can be provided by 
simple analytical means (Girimaji, 1991b). 

Finally, in the context of single-point PDF formulation presented, there is no information 
pertaining to the evolution of the relevant turbulent length scales. The final expressions can 
be only presented in terms of other physical parameters (here, only through the parameters 
g and/or £). In the context considered, this parameter has been provided by the DNS data. 
In a practical application, this must be provided by external means (turbulence models, 
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experimental data, etc.) Jiang et al. (1992), and Frankel et al. (1992) provide further 
discussions related to this issue. 
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5 Concluding Remarks 

It is demonstrated that the mapping closure of Kraichnan-Pope (Chen et al ., 1989; Pope, 
1991) yields closed form analytical expressions for predicting the limiting bound of reactant 
conversion rate in a simple chemistry of the type F + rO — *• (1+r) Product in homogeneous, 
isothermal turbulence. It is also shown that for the case of complete initial segregation, 
the scalar PDF’s generated by this closure can be well-approximated by a Beta density. 
This density also provides closed form analytical expressions for the limiting rate of reactant 
conversion. A nice feature of the mathematical results generated by the two models is their 
capability of revealing the influence of the stoichiometric coefficient and the equivalence 
ratio. In both cases, the mathematical expressions simplify significantly for a stoichiometric 
mixture. The prediction results via both models compare favorably with data generated by 
DNS. This agreement is not surprising, since for the case of a mixing controlled reaction, it 
follows from the compatibility of the models PDF’s with those of DNS. However, the simple 
final results generated here are superior to those of previous closures based on the typical 
C/D models, and those based on Toor’s hypothesis. For example, in comparison with Toor’s 
results, it is shown that in stoichiometric mixtures, an improvement of approximately 36% 
can be made in evaluating the unmixedness; for non-stoichiometric mixtures the improvement 
can be as high of 100%, depending on the value of the equivalence ratio. 

These closed form relations are furnished with the imposition of several restrictive assump- 
tions. The ramifications associated with these assumptions are discussed, and some sugges- 
tions for future extensions are provided. In spite of these assumptions, it is very encouraging 
to have physically plausible algebraic relations for direct estimate of the reactant conversion 
rate in plug flow reactors. Because of their simplicity and demonstrated validity, these ex- 
pressions are strongly recommended for routine and economical engineering predictions in 
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unpremixed binary reacting systems such as those in plug flow reactors. 


Nomenclature 

a, 6, c. some constant. 

B. the Beta function. 

V. the Parabolic Cylinder Function. 

Da. Damkohler number. 

F. Fuel. 

Q. parameter in the mapping closure. 

g. half the inverse of the normalized variance of the Shvab-Zeldovich variable. 
H. the Heaviside function. 

T. the incomplete Beta function. 

J . Shvab-Zeldovich variable. 

0 . Oxidizer. 

V. PDF. 

r. the stoichiometric coefficient. 
t. physical time 

W. area weight of the reactant. 
y. dummy variable of integration. 

Z. the unmixedness ratio. 

Greek letters 

/?!,/? 2, parameters of the Beta density. 

7. equivalence ratio. 
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T. the Gamma Function 
6. the delta function. 

the composition space for the PDF. 

<f> 0 . composition space for a Gaussian reference field, 
r. normalized time. 

X- the mapping function. 

Subscripts 

0. time zero (inlet of plug flow reactor). 

G. Gaussian, 
st. stoichiometric. 


Other symbols 

( ). Probability Average. 
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Figure Captions 


Figure 1: PDF of the Shvab-Zeldovich variable at two times ( t2 > <1) for the symmetric case 

(<j>-j). 

Figure 2: PDF of the Shvab-Zeldovich variable at two times (t2 > tl) for the non-symmetric 
case (< J >= 0.625). 

Figure 3: Normalized mean concentration of fuel and the oxidizer for the symmetric case, 
and Ja t = 0.5. 

Figure 4: Normalized mean concentration of fuel and the oxidizer for the non-symmetric 
case, (a) J at = 0.4, (b) J at = 0.2. 

Figure 5: The unmixedness ratio at several values of the equivalence ratio. 
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Results are presented of direct numerical simulations (DNS) of a spatially developing shear flow 
under the influence of an infinitely fast chemical reactions of the type A + B — *■ Products. The 
simulation results are used to construct the compositional structure of the scalar field in a statistical 
manner. The results of this statistical analysis indicate that the use of a Beta density for the probability 
density function (PDF) of an appropriate Shvab-Zeldovich mixture fraction provides a very good 
estimate of the limiting bounds of the reactant conversion rate within the shear layer. This provides a 
strong justification for the implementation of this density in practical modeling of non-homogeneous 
turbulent reacting flows. However, the validity of the model cannot be generalized for predictions of 
higher order statistical quantities. A closed form analytical expression is presented for predicting the 
maximum rate of reactant conversion in non-homogeneous reacting turbulence. 


KEYWORDS Direct numerical simulations Turbulent shear flows Beta density 
Probability density function Unmixcduess. 


INTRODUCTION 

In our previous work (Madnia et al , 1991a) we have demonstrated that the Beta 
density provides a very good means of approximating the probability distribution 
of a conserved scalar quantity in homogeneous turbulent flows. The versatility of 
this density enables it to predict the evolution of a conserved scalar variable in a 
statistically homogeneous field, in a manner that is in accord with both DNS and 
experimental measurements. Despite this positive feature, some questions 
pertaining to the generality of this density for more realistic flow configurations 
remained unanswered. 

In this work, we intend to address some of these remaining questions. 
Specifically, we shall make a comparative assessment of the behavior of this 
model, and provide an estimate of the degree of its validity in the context of more 
general conditions than previously considered. For this purpose, we have selected 
a configuration of particular interest in the design of diffusion controlled reactors; 
namely, a spatially developing parallel mixing layer under the influence of 
harmonic perturbations. The flow field which develops in this setting is dominated 
by large scale coherent structures and exhibits non-homogeneous features. This 
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inhomogeneity plays a dominant role in the mechanism of mixing in turbulent 
flows and is very suitable for promoting mixing and chemical reactions. Because 
of this property, this configuration has been widely utilized in numerous 
investigations on turbulent reacting flow phenomena (see Givi and Riley (1992) 
for a recent review). 

We shall follow an analogous approach to that of our previous work (Madnia et 
al.y 1991a) in assessing the model's behavior. However, there are fundamental 
differences between the works in formulating the problem and in implementing 
the numerical procedure involved in the DNS. The configuration considered here 
is non-homogeneous. Therefore, both the procedure of statistical sampling and 
the numerical algorithm used in the DNS are different from the homogeneous 
simulations of Madnia et al. Also, the mathematical operations leading to the 
estimate of the reactant conversion rate are substantially more elaborate. 
However, analytical solutions can be attained with the aid of the mathematical 
properties of the Beta Function and the Incomplete Beta Function. With these, 
the generated results can be directly compared with those of DNS. 


DESCRIPTION 


The problem under consideration is that of a spatially evolving mixing layer 
containing initially segregated reacting species A and B in the two free streams. 
Species A is introduced into the upper, high-speed stream and species B enters on 
the lower, low-speed side. The chemical reaction between the two species is 
assumed to be single step, irreversible and infinitely fast. The magnitude of the 
Mach number is assumed negligible. Therefore, the flow is considered incompres- 
sible. The two species are assumed thermodynamically similar with identical 
diffusion coefficients. Within the framework of these approximations, the flow is 
mathematically described by the conservation equations of mass, momentum, and 
a species conservation equation for the trace of a conserved Shvab-Zeldovich 
variable /, which characterizes the compositional structure of the flow. These 
equations are parameterized by the Reynolds number, the Peclet number, and 
the velocity ratio across the layer. 

The computational package employed in the DNS is based on a hybrid 
pseudospectral-spectral element algorithm developed by Givi and Jou (1988) and 
McMurtry and Givi (1991). The pseudospectral routine, which is based on 
Fourier collocation, is used in the cross stream direction together with free slip 
boundary conditions. The spectral-element discretization is employed in the 
streamwise direction. This involves the decomposition of the domain in this 
direction into an ensemble of macro finite elements. Within each of these 
elements, the dependent flow variables arc approximated spectrally by means of 
Tchebysheff polynomials of the first kind. The hybrid procedure employed in this 
way is very attractive in that it combines the accuracy of spectral discretization 
with the versatility of finite element methods. Therefore, it is convenient for 
accurate simulations of the complex spatially developing flow under consideration 
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here. For a detailed description of this hybrid method for DNS of turbulent shear 
flows, the reader is referred to McMurtry and Givi (1991). 

The flow field is initialized at the inflow by a hyperbolic tangent mean velocity 
distribution, upon which low amplitude perturbations are superimposed in order 
to promote the formation of coherent vortices. The initialization of the 
Shvab-Zeldovich variable at the inflow is such that it takes on the values of 
zero and one in the streams containing B and A, respectively. The disturbances 
on the mean flow correspond to the most unstable mode of the hyperbolic tangent 
profile (Michalke, 1965) and its first two subharmonics. The amplitude of the 
fluctuating velocity is set equal to 6% of the mean velocity. This amplitude 
corresponds to that measured in typical shearlayer experiments. The magnitudes 
of the phase shifts between the modes of the instability waves are randomly 
selected from a random seed with a top hat PDF of zero mean and a specified 
variance. The implementation of these phase shifts is the only mechanism to 
introduce randomness into an otherwise deterministic simulation. This is to 
partially mimic the random “commotions” which are present in the “universe,” 
into an isolated two-dimensional simulation. At the outflow, a weak condition of 
zero second derivative is applied for all the dependent variables. This boundary 
condition cannot be justified in a rigorous mathematical sense and can be 
specified only in an ad hoc manner. This was done in order to facilitate 
implementation of the Dirichlet boundary conditions in the finite element 
procedure employed here. Moreover, the results of extensive numerical tests 
showed that the effects of the approximate boundary conditions are confined 
within the last computational element. This is consistent with that found 
previously by Korczak and Hu (1987). Therefore, the solution in the last 
computational domain can be ignored. 

With the solution of the unsteady transport equations, the DNS data are 
extremely useful in visualizing the instantaneous behavior of the flow. The results 
of these direct simulations can also provide an assessment of the statistical 
behavior of the flow. This is due to the availability of flow information at all the 
computational grid points and at all times. Therefore, statistical sampling can be 
done quite easily by ensemble averaging of the instantaneous data gathered over 
many realizations. 

In our simulations the instantaneous data of the Shvab— Zeldovich variable, 
obtained from the DNS, are used for statistical analysis. A total of 3200 
realizations were employed in the sampling. With this ensemble, the magnitudes 
of the mean and the variance of the Shvab— Zeldovich variable are calculated at 
all the grid points. These statistical quantities are of primary interest for 
turbulence modeling. Based on the knowledge of these first two moments, a Beta 
density is assumed for the PDF of the Shvab-Zeldovich variable. The ensemble 
mean values of the product concentration and the unmixedness are calculated 
subsequently via this density. These modeled quantities are then compared with 
those obtained directly from the DNS data. At the majority of grid points the 
Beta density is asymmetric. Nevertheless, analytical solutions for the statistics of 
the reacting field are possible and have recently been obtained by Madnia et al. 
(1991b). The derivations of the equations are rather involved. Here, we only 
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present the final results for the normalized mean product concentration (C p ) and 
that of the unmixedness X V 2 : 


< C p ) = (C p )(x,y) = l-e l -8 2 ( 1 ) 

'P 2 = 'P 2 (x,y) = -0 1 d 2 (2) 


where: 


R 1 r(o- + ff) a-ff . , R ., 

1 2 a+t, ~\a+P)r(a)r(p) a + /3 (1 l **' ' P)) 


R 1 T{a+P) B-a 

2 2 a+p ~\a + p) r(or)r08) a + P *" K ’ P) 


(3) 

(4) 


Here, T is the Gamma function, and oc and /? are related to the first two moments 
of the random field (Frankel et al y 1991). £) is the Incomplete Beta 

Function, and J> sx denotes the stoichiometric value of the Shvab-Zeldovich 
variable, For unity normalized concentrations at the free streams, / st = 0.5. 
The results predicted by these equations shall be compared with those generated 
by DNS for a quantitative assessment of the closure. 


PRESENTATATION OF RESULTS 

Computations were performed in a domain of size 135 x 346, where 6 denotes 
the vorticity thickness at the inlet of the flow. There are 65 Fourier modes in the 
cross-stream direction and 42 finite elements arc used for the strcamwisc 
discretization. A fifth order Tchebysheff polynomial is employed to approximate 
the variables within each element. This discretization is equivalent to, at least, a 
fifth order accurate finite difference technique even if the spectral convergence is 
not considered. This results in a total number of 211 points in the streamwise 
direction. A second-order Adams-Bashforth finite difference scheme was utilized 
for temporal discretization. With the computationally affordable 211x65 grid 
resolution available, accurate simulations with Re = Pe = 200 (where these 
normalized quantities are defined based on the velocity difference across the 
layer, and the initial vorticity thickness) are possible (Givi, 1989). These values 
are somewhat smaller than those of a fully developed laboratory turbulent shear 
layer, but are sufficiently high to promote the growth of instability waves. 

In order to visualize the flow evolution, a time series of DNS generated 
snapshots of the instantaneously normalized product concentration are shown in 
Figure 1. This figure depicts the way in which the small amplitude perturbations 
manifest themselves in the formation of large scale coherent vortices downstream 
of the splitter plate. The forcing associated with the most unstable mode alone 
would cause the initial roll-up of these vortices. The perturbations corresponding 
to the first subharmonic mode result in a second roll-up in the form of merging 
neighboring vortices, and the presence of the second subharmonic generates a 
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FIGURE 1 Time sequence of contour plots of instantaneous product concentration. 

(See color Plate I.) 

third roll-up (second merging) at a region near the outflow. The roll-up and the 
subsequent pairings of these vortical structures engulf the unreacted species from 
the two streams into the mixing zone. The dynamics of the vortical evolution in 
this process reflect the two main mechanisms of mixing enhancement. Large 
concentrations of vorticity, that are approximately of one sign, bring the fluid 
from the two free streams into the large scale structures. Within these structures, 
the fluid elements are subject to further straining, and diffusion between the two 
fluids results in final mixing at the molecular scale. With these processes, along 
with the assumption of infinitely fast chemistry, the rate of product formation is 
naturally enhanced by the dynamics of vortical evolution. 

The ensemble average of the product concentration, formed as a result of the 
mechanisms just described, is presented in Figure 2(a). The statistical averaging 
process was performed over a time span approximately equal to the residence 
time of the flow within the domain of computational consideration. This analysis 
was performed after the effects of the initial transients were washed away at the 
outflow. A comparison between Figures 1 and 2(a) reveals the essential 
deficiencies of tile ensemble averaged results. Note that, all the detailed 
intricacies of the dynamics of flow evolution are masked by the averaging. While 
these ensemble-averaged results are of main interest for practical applications. 
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FIGURE 2 Contour plots of mean product concentration generated by (a) DNS, (b) the model 

(See color Plate II.) 


they are not capable of describing many of the interesting features of turbulence 
transport. The results depicted in Figure 2(a) are the best that one can expect to 
generate from a turbulence closure in a statistical sense. 

To demonstrate the capability of the model based on the Beta density, in 
Figure 2(b), the average product concentration obtained by means of Eqs. 
(l)-(4) is presented. A comparison of this figure with Figure 2(a) reveals a 
remarkable similarity between the two results. For a more quantitative com- 
parison, the cross stream variation of this product concentration and those of the 
negative of the unmixedness are presented at two streamwise locations in Figures 
3-4. Also the streamwise variation of the total product, T p , is presented in Figure 
5. In all these figures, the comparisons between the model predictions and the 
DNS results are noteworthy. 

The agreements demonstrated by the comparisons noted above provide a 
reasonable justification for recommending Eqs. (l)-(4) as working relations in 
routine engineering predictions. However, these relations are advocated only for 
statistical predictions of low order moments. The physics of the mixing phenome- 
non in a spatially developing shear flow is far too complex to be completely 
described by the Beta density. The highly intermittent nature of the mixing 
process exhibits features that cannot be fully captured by this density (Givi, 
1989). To demonstrate this point quantitatively, in Figures 6-7 the cross stream 
variations of the third and fourth centralized moments of the Shvab Zeldovich 
variable calculated from the Bela density are compared with those generated by 
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Collocation points 


FIGURE 5 Streamwise variation of total product (7^) from the DNS and the model. 

DNS. These figures indicate that while the trend is somewhat similar in the two 
cases, the extent of agreement is not as good as those of lower level statistical 
quantities (Figures 2-5). This is primarily due to the mechanism of large scale 
entrainment in the shear flow, in that there are always unmixed fluids present 
within the core of the vortical structures. This results in a slight trimodal behavior 
in the probability distributions generated by DNS (Givi and Jou, 1988; Lowery 
and Reynolds, 1990) and observed experimentally (Koochesfahani, 1984, Dimota- 
kis, 1989). This trimodal behavior cannot be captured by the Beta density which 
is at best suitable for producing bimodal distributions. 

Recall that the specification of the Beta density is based on the knowledge of 
the first two moments of the conserved scalar. Here, these moments were 
supplied via DNS, but in an actual engineering application they must be provided 
by other means ( i.e . an appropriate turbulence closure, or experimental data). 
However, specification of these parameters is substantially simpler than those of 
the reacting scalars, even with the assumption of an infinitely fast reaction. This is 
simply due to the fact that such data can be obtained in simpler, non-reacting flow 
experiments without any physical complications associated with the chemistry. 
With such data available, then Eqs. (1)— (4) are deemed very suitable, at least in 
the absence of better alternatives, for a reasonably accurate and inexpensive 
prediction of the limiting bounds of the reactant conversion rate. These equations 
are therefore recommended for this purpose in either homogeneous or non- 
homogeneous flows. 
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NOMENCLATURE 

A y B Reactants. 

C p Product Concentration. 

$ Shvab-Zeldovich variable. 

Pe The Peclet number. 

Re The Reynolds number. 

T p (x) Total product = j-Zdy j x 0 dx{C p )( x, y). 
x ,y The physical coordinates. 

<5 Vorticity thickness at the inflow. 

W 2 The unmixedness, 

ju 3 Third order moments of the Shvab-Zeldovich variable. 

Fourth order moments of the Shvab-Zeldovich variable. 

() Ensemble average. 
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Results are presented of direct numerical simulations (DNS) of a randomly per- 
turbed compressible, spatially developing planar jet under the influence of a finite rate 
Arrhenius chemical reaction of the type F + O — * Product + Heat. The objectives 
of the simulations are to assess the compositional structure of the flame, and to de- 
termine the influence of reaction exothermicity by means of statistical sampling of 
the data generated by DNS. It is shown that even with this idealized kinetics model, 
the simulated results exhibit features in accord with experimental data. These results 
indicate that the Damkohler number is an important parameter in determining the 
statistical composition of the reacting field, and that the results are insensitive to the 
mechanism by which this parameter is varied. It is demonstrated that as the intensity 
of mixing is increased and the effect of finite rate chemistry is more pronounced, the 
magnitudes of the ensemble mean and rms of the product mass fraction decrease, and 
those of the reactants mass fraction Increase. Also, at higher mixing rates the joint 
probability density functions of the reactants’ mass fractions shift towards higher val- 
ues within the composition domain indicating a lower reactedness. These trends are 
consistent with those observed experimentally and are very useful in portraying the 
statistical structure of non-equilibrium diffusion flames. The DNS generated data are 
also utilized to examine the applicability of the “laminar diffusion flamelet model” in 
predicting the rate of the reactant conversion with finite rate chemistry. This exami- 
nation indicates that the bounds of the product formation scale reasonably well with 
those obtained by the flamelet model. However, in the range of the relatively low values 
of the Damkohler numbers considered, the scatter of the results is substantially more 
than that to be completely determined by the closure. Finally, the simulated results 
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suggest that in the setting of a “turbulent” flame, the effect of the heat liberated by 
the chemical reaction is to increase the rate of reactant conversion. This finding is 
different from that of earlier DNS results and laboratory investigations which indicate 
a suppressed chemical reaction with increasing heat release. 
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1 Introduction 

Recent advances in flow visualization and diagnostic techniques have made significant con- 
tributions to the understanding of turbulent combustion phenomena [1, 2]. With implemen- 
tation of these advanced techniques it is now possible to make a detailed examination of 
the structure of turbulent flames, and to establish a better understanding of the intricate 
physics of such flames. Amongst recent applications of these diagnostic techniques, Masri et 
al. [3, 4] have made noteworthy progress in describing the structure of unpremixed turbulent 
jet flames under the influence of finite rate chemical reactions (see [5, 6] for comprehensive 
reviews). In particular, they have reported the results provided by statistical sampling of 
the data obtained by Raman- Rayleigh scattering measurements of the species mass fraction 
and the temperature in an unpremixed turbulent methane jet flame. These measurements 
are utilized primarily for constructing single-point statistics of the reacting field, includ- 
ing the ensemble mean, rms and the probability density functions (PDF’s) of the relevant 
thermochemical variables. The response of the flow under different mixing conditions is 
characterized by a comparison of the results obtained at various mixing rates. The sampling 
of the data in this way allows a systematic assessment of the compositional structure of the 
non-equilibrium flame as it approaches extinction. Such data have proven very valuable in 
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describing the phenomena of mixing and chemical reactions in unpremixed turbulent com- 
bustion, and in addressing some of the existing problems in the modeling of such flames by 
means of statistical methods. 

In order to provide a computational complement to laboratory investigations, there has been 
significant progress in the development and implementation of direct numerical simulations 
(DNS) of unpremixed reacting flows [7, 8, 9, 10, 11, 12, 13, 14] (see [15, 16, 17, 18, 19] 
for recent reviews). The extent of progress made by these contributions has been very 
encouraging, thus justifying further utilization of direct methods in the analysis of turbulent 
flames. In this work, we intend to continue our investigation of turbulent flames via DNS, 
and to simulate such flames in a context relevant to laboratory investigations. Obviously, it 
is still impossible to perform DNS of a realistic “physical” flame [15]. Here, our intention 
is to make use of the present capabilities of supercomputers to further address some of the 
fundamental and pertinent issues in regard to the dynamics of such flames. 

In the next section, the method of investigation is presented together with the underlying 
assumptions imposed in making the problem tractable for computational treatment. The 
ramifications associated with the simplifying assumptions are discussed, where appropriate, 
in the presentation of results. This presentation is made in Section 3 with a focus on 
assessing the compositional structure of the flame, and on revealing the influences of reaction 
exothermicity on the evolution of the flame. With this presentation, it is shown that some 
of the trends observed in laboratory experiments can be captured by DNS, and speculation 
is made on some of the other features not yet considered in laboratory investigations. This 
paper is drawn to a conclusion in Section 4. 
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2 Method of Investigation 


The flow under investigation consists of a spatially developing, two-dimensional, planar jet 
flame. The schematic diagram of the jet is shown in Fig. 1(a) together with the specifications 
of the physical dimensions. The fuel, F , discharges from the inner higher speed jet of width 
D with a velocity Uf into a co-flowing lower velocity ( Uo ) oxidizer, O, stream. The inlet 
temperature and pressure are the same at the inlet of both streams, and the reactants are 
completely segregated at the inlet. This configuration is different from that of a circular 
jet flame considered in [3, 4], but allows the elucidation of some specific physical features 
that are believed to be invariant of the jet geometry [18]. The flow evolves spatially in the 
streamwise direction (x), and impermeable boundary conditions are imposed in the cross- 
stream direction ( y ). The DNS computational procedure is the same as that employed in 
[14, 20]. This involves the discretization of the complete Navier-Stokes equations as well 
as the applicable Fickian chemical species conservation equations. This discretization is 
based on a two-four (second order accurate in time, and fourth-order truncation in space) 
compact parameter finite difference scheme. This scheme has been developed and made fully 
operational by Carpenter [21] and is utilized in the context of the SPARK computer code 
developed by Drummond [22]. A generalized single step Arrhenius exothermic chemistry 
model of the form F + O —* Product + Heat is assumed. This kinetic mechanism is simple 
enough to allow an efficient treatment via DNS, and yet it is capable of capturing some of the 
important non-equilibrium effects associated with combustion. All the species involved in 
this reaction are assumed to be thermodynamically identical, and the values of the physical 
transport parameters are assumed to be invariant with temperature. With this kinetics 
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prototype, the rate of reactant conversion is expressed by: 



u = K f C F C 0 exp(-—) 


( 1 ) 


iu The flow field is initialized with a top hat streamwise velocity profile with a low amplitude 

1 random forcing at the jet entrance. The magnitudes of the forcing frequency and the phase 

L 

shifts between the different components of the forcing eigenfunctions are selected randomly 
*■" from a seed with a uniform probability distribution and with specified values of the mean 

and variance. This random specification of the frequency and the phase shift is intended to 
mimic the influence of turbulence in an otherwise deterministic simulation. The amplitudes 
of the perturbation are chosen so that the maximum value of the turbulent intensity at the 
— inflow is no more than 5%. 
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The parameters influencing the rate of reactant conversion in this setting are the Damkohler 
number, the Zeldovich number, the Reynolds number, the heat release parameter, the con- 
vective Mach number, the velocity ratio between the streams, and the Prandtl/Schmidt 
numbers (see the nomenclature for a description of the thermochemical variables and for a 
definition of all the relevant non-dimensionalized parameters). Based on our earlier find- 
ings [20, 14], in the range of parameters considered the Reynolds number does not have a 
significant effect. Therefore, the influences of this parameter as well as those of the Peclet 
and the Schmidt numbers are not investigated. The values of the molecular Prandtl and 
the Schmidt numbers are set to unity ( Pr = Sc = 1), to make the analysis associated 
with laminar diffusion flamelet modeling easier. The magnitude of the Zeldovich number 
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is held constant in all the simulations, but the values of the Damkohler numbers and the 
heat release parameter are altered to assess their influence on the dynamics of combustion. 
The full compressible form of the transport equations are considered in the computational 
formulation. However, the magnitudes of the convective Mach numbers considered are not 
very large. Therefore, many of the physical issues associated with high speed combustion 
[22, 18, 14] are not addressed. Also, the body forces are assumed negligible; therefore, the 
interesting physical problems associated with the buoyancy [23, 24] are not considered. 

The computer code used is ideal for high resolution, three-dimensional (3-D) simulations in 
a reasonably efficient manner on presently available supercomputers [16, 20]. However, the 
specific physical issues considered here do not require 3-D treatment, and can be addressed 
in a 2-D context. This was established by some comparative assessments of 2-D and 3-D 
results, and all subsequent simulations were performed in a 2-D domain. This yields a higher 
numerical resolution at a lower computational cost. 
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In accordance with the procedure implemented in [3], three flames are considered; Flames 
A, B, and C. These flames are identified by the relative magnitude of the mixing intensity 
quantified by the value of the Damkohler number. In Flames B and C, the magnitudes of 
the Damkohler number are the same and are set equal to half of that in Flame A. However, 
the mechanism by which this is enacted is different. In Flame B, this is implemented by 
decreasing the magnitude of the chemical reaction frequency, and in Flame C by increasing 
the magnitude of the characteristic velocity of the jet. Both of these procedures result in the 
same reduction in the value of the Damkohler number, and the difference between Flames A 
and C is the same as that in the flames considered in [3]. The magnitude of the convective 
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Mach number is higher in Flame C than that in Flame A. But, the value of this Mach 
number is not so large as to have a substantial influence on the spatial growth of the jet [14]. 
Also, the imposition of external forcing masks the effects of other possible mixing inhibition 
mechanisms (such as exothermicity, density stratification, etc.) [12, 24, 14]. 

The compositional structure of the flame is assessed by means of examining the response 
of the flow to the Damkohler number. This assessment is quantified by considering the 
statistical properties of the reacting field, including the ensemble mean, rms and the PDF’s of 
the mass fractions of all the species. The analysis has been made first in a non-heat releasing 
reacting jet. The important influences of the exothermicity are subsequently considered. 


3 Presentation of Results 

With available computational resources, numerical experiments with a resolution of 245 (in 
a;) xl65 (in y) mesh are possible. This resolution is sufficient provided that an appropriate 
grid generation procedure with a non-uniform mesh and with a heavy concentration of grid 
points near the regions of maximum instantaneous shear is utilized (Fig. 1(b)) [22, 20, 14]. 
With this configuration, it is possible to perform simulations with moderate values of the 
Reynolds, Peclet, Damkohler and Zeldovich numbers within a domain L x = 13 \D > x > 
0, L y = 6§D > y > 0. With the numerical resolution attained, it was possible to perform 
simulations with a Reynolds number of Re 6 x 10 4 . The magnitude of the Damkohler 
number is Da = 5 for Flame A and Da = 2.5 for Flames B and C. The Zeldovich number is 
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set at Ze = 20 for all the flames, and the heat release parameter is in the range Q = 0 ~ 5. 
The velocity ratio between the two streams is With the resolution adopted, a typical 

simulation requires about 6 hours of CPU time on a CRAY-2 supercomputer. 

To visualize the global structures and the unsteady evolution of the layer, a time sequence 
of the DNS generated instantaneous product mass fractions is shown in Fig. 2(a). This 
figure shows how the growth of low amplitude random perturbations manifest themselves in 
the formation of large scale coherent vortices downstream of the jet entrance. Due to the 
random nature of the forcing mechanism, the evolution of the large scale structures is very 
complex and, similar to laboratory shear flows, consists of random formation of vorticity 
rollups and subsequent pairing and coalescence of the neighboring vortices. Also, due to 
this randomness, the layer is asymmetric with respect to the jet centerline even though a 
symmetric disturbance is imposed at the inflow. In all the simulations, the energy spectrum 
of the fluctuating hydrodynamics field is fairly broad-banded. Therefore, the flow can be 
assumed “turbulent”, albeit two-dimensional. The formation of large scale structures results 
in an enhancement of mixing, as demonstrated in Fig. 2(a). The unsteady conformation of 
large scale structures results in the mutual engulfment and entrainment of the two reactants, 
and the dynamics of this evolution reflects the primary mechanisms of the mixing process. 
Large concentrations of vorticity on both sides of the jet centerline, bring the fluid from 
the two jets into the core of the large scale structures. Within these structures, the fluid 
elements are subject to further straining, and the molecular diffusion between the two fluids 
results in final mixing at the molecular scale. The enhanced mixing mechanism caused by 
the large scale vortices results in a general augmentation of the product formation. At 
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regions close to the jet entrance, the chemical reaction is promoted by the formation of 
coherent vortices, and there is a substantial amount of product formed. If the flow was in 
chemical equilibrium, the total amount of product would have increased monotonically along 
the streamwise direction. For the non-equilibrium chemistry considered here, this is not the 
case and mixing enhancement does not necessarily yield improved combustion. It is shown 
in Fig. 2(a) that at the shear zone near the entrance of the jet, where the reactants are first 
brought into contact, the extent of chemical reaction is fairly uniform. Further downstream, 
as the magnitude of the instantaneous strain increases, the reaction rate approaches zero at 
the braids of the vortices and the flame locally quenches at those locations. This mechanism 
of flame extinction is consistent with previous DNS results [9, 10], and is also corroborated 
by experimental observations [25, 26]. This behavior was observed in all the simulations 
with moderate Damkohler numbers asserting the departure from chemical equilibrium. 

The effects of the non-equilibrium chemistry become more clear as the influence of finite 
rate kinetics is more pronounced. To show this influence in the format adopted in [3], the 
results extracted by DNS are statistically analyzed. With the assumption of statistically 
stationary flow, the instantaneous data can be used to construct ensemble-averages. Here, 
these averages are obtained by statistical analysis of 1250 time samples gathered after the 
effects of initial transients are washed away at the outflow. The data sampling was performed 
within a time period equivalent to half the residence time of the flow, which was shown to 
yield statistically acceptable trends. These results are presented in Fig. 2(b) in the form 
of contour plots of ensemble- averaged product mass fraction. A comparison between the 
two parts of Fig. 2 reveals the essential deficiencies of the ensemble averaged results. Note 
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that the intricate dynamics of the flow evolution are masked by the averaging procedure. 
While these ensemble- averaged results are of interest for practical applications, they are not 
capable of describing the interesting detailed features of turbulent transport. Nevertheless, 
the results presented in Fig. 2(b) are those which one can expect to generate by using a 
turbulence closure. Therefore, an assessment of their validity is essential for the purpose of 
statistical modeling. These results also facilitate the most convenient (and perhaps the most 
traditional) means of comparing our DNS data with experimental measurements. 

To exhibit the effects of finite rate chemistry in a quantitative manner, the statistical analysis 
of the data is very useful. With DNS data available at every time step at all the locations 
within the flowfield, an extensive library of data sets is produced. In the discussions be- 
low, the statistical data are presented at several streamwise locations where the appropriate 
physics are captured. In Fig. 3 the cross stream variation of the ensemble-mean values of the 
product mass fraction is presented at a streamwise location (x = 10D). This figure shows 
that the magnitude of mean product mass fraction in Flame A is substantially higher than 
those of Flames B and C. This suggests that as the flame approaches extinction, the mass 
fractions of the reactants in the reaction zone increase and, as a result of local flame quench- 
ing, the rate of product formation decreases. The same behavior is observed in Fig. 4, which 
shows that the amplitude of the rms of the product mass fraction is higher in Flame A than 
those in Flames B and C. The trends observed in these two figures are consistent with those 
of the laboratory measurements [3] and reveal the effects of the mixing intensity on com- 
positional structure of the flame. These results also indicate that this behavior is observed 
regardless of the procedure by which the Damkohler number is varied. In the experiments, 
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the Damkohler number was altered only by means of increasing the hydrodynamic intensity 
(be. increasing the inner jet velocity). This was because a single fuel was used and, there- 
fore, the chemical frequency was constant in all the experiments. Our results indicate that 
this dependence is observed regardless of the procedure by which the change is made, and 
Flames B and C portray approximately the same statistical conduct. The behavior observed 
in Figs. 5-6 demonstrates the dominant influence of the Damkohler number in quantitative 
analysis of non-equilibrium flames, and also in the comparative assessment of the statistical 
results generated in our DNS with those obtained in laboratory experiments [3] (See [27] for 
more discussions pertinent to this issue). 


i m 




The influence of finite rate chemical reactions on the compositional structure of the flame is 
further demonstrated by examining the marginal and the joint PDF’s of the reacting scalar 
variables, as well as the “scatter” plots of the product mass fraction. These are presented 
in Figs. 5-8. In Fig. 5, results are presented of typical joint PDF’s of the Fuel and the 
Oxidizer mass fractions, be. V{Yp,Yo) for Flames A and B. This figure indicates that at 
lower mixing rate, Flame A (Fig. 5(a)), there is a larger amount of product formed, and 
the height of the PDF’s are larger at the mixed mass fraction values at the interior of the 
allowable region (0 < Y A ,Y B < 1). At the higher mixing rate, Flame B (Fig. 5(b)), the 
mixture is composed of the two reactants with the joint PDF’s concentrated primarily along 
the line of Yp + Yo = 1. These results are consistent with those reported in [4] and indicate 
that at the higher levels of mixing intensity, the results are closer to those obtained by cold 
mixing. 


The corresponding marginal PDF’s of the product mass fraction V(Y p ) are shown in Figs. 6 
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for Flames A and B. In both flames, the PDF’s are concentrated near zero at free streams, 
and become broader near the central region of the layer. As the mixing rate is decreased, the 
PDF’s are shifted toward higher values of the mass fraction. This behavior is also consistent 
with the experimental data on PDF’s of product species, H ? , CO 21 • • • [4] an d indicate lower 
reactedness at higher mixing rates. However, due to strong intermittency of the flow in 
the two-dimensional simulations, the PDF’s do not adopt a Gaussian distribution, and are 
always heavily populated near free stream values, i.e. zero. Also, the tails of the PDF’s get 
thinner as the mixing is intensified. This suggests a reduction in the rms of the product mass 
fraction, as was reported in the experiments of Masri et al. [3], and also evident in Fig. 4. 

The lesser rate of product formation is also evident in the scatter plots of the product mass 
fraction vs. the mixture fraction in Figs. 7 and 8. These figures represent conditions at 
several streamwise locations, including all the values in the cross stream direction. There 
are a total of 102500 data points taken in devising these scatter plots. It is indicated in these 
figures that the rate of product formation is generally high near the inflow, and gradually 
decreases further downstream as the influence of flame extinction becomes more pronounced. 
Note that at x = 6.8D, the reaction zone is confined to a small region near the equilibrium 
limit, but at distances further downstream the flame becomes widely distributed within the 
allowable region. This spread is more noticeable for the lower Damkoher number flame, 
indicative of the pronounced effects of the extinction. This is to be expected, since in the 
context of the simplified kinetics model, it is expected to have an infinitely thin reaction zone 
near the equilibrium limit, as the magnitude of the Damkohler number becomes infinitely 
large. With a departure from this limiting rate, the scatter becomes wider as the Damkohler 
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number is reduced. Also, as the influence of flame extinction becomes more pronounced (at 
lower Damkohler numbers), it is not surprising to have a lower rate of product formation, 
as evidenced by the smaller amplitudes of the mass fraction values in Fig. 8 in comparison 
with those in Fig. 7. 

In order to have a quantitative means of departure from the equilibrium, the simulated results 
in Figs. 7-8 are also compared with those predicted by the steady laminar diffusion flamelet 
model proposed in [26, 28]. For this comparison, the simplified chemical kinetics model 
employed above was also considered in the content of a simple one-dimensional opposed 
jet laminar flame configuration similar to those in [29, 30, 31]. In this setting, with the 
assumption of unity Lewis number, the simulations based on the flamelet model involves the 
solution of the diffusion-reaction equation of type [31]: 

(PjP __ ~ 2u y 
RX ' 

In this equation, 0 denotes the reacting scalar variable, i.e. species mass fractions, and the 
temperature, with the corresponding chemical source term, uty. The variable \ represents 
the dissipation of the conserved mixture fraction, £, is expressed by: 


X = Xst exp 


-2 (erf -1 (2£ — 1))' 


( 3 ) 


Here, Xst denotes the “stagnation” value of the scalar dissipation. In the laminar jet flame, 
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this stagnation value is determined directly by the hydrodynamic intensity of the jet. How- 
ever, for the unsteady, “turbulent” jet flame considered here, a single value cannot be 
adopted. Instead, only some representative values are considered. The results correspond- 
ing to several of these values are superimposed on Fig. 7. Obviously, as the magnitude of 
the stagnation dissipation increases, the rate of fuel consumption decreases, and the results 
are always bounded in the upper limit as the value of the Damkohler number approaches 
infinity. Now, holding the value of the stagnation dissipation constant together with the 
other important kinetics parameters, but decreasing the value of the Damkohler number, 
the results can be superimposed on the scatter plots in Fig. 8. With the reduction of the 
Damkohler number the rate of product formation becomes less, and the laminar diffusion 
flamelet model shows similar trends as those of the DNS data. However, as shown in both 
Figs. 7 and 8, with a fixed pre-determined value of the stagnation dissipation it is not possible 
to capture the complete picture, as the predicted results cannot account for the wide scatter 
in the figures. The comparison between the model predictions and the DNS results would 
probably be better at higher Damkohler numbers, with a thinner reaction zone, and perhaps 
by employing a “variable” form for the stagnation scalar dissipation. It is not presently clear 
how to modify this rate to account for the strong effects of the spatial inhomogeneities, and 
it is speculated that the functional form of the stagnation dissipation is strongly dependent 
on the details of the flow field. Nevertheless, the simulated results with this strong finite rate 
kinetics effects are consistent with the laboratory data, as similar behavior was also observed 
in the experimental measurements [4]. 

Despite the attractive features of our simulations in predicting the trends observed experi- 
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mentally, it must be emphasized that in the context considered, the factors influencing the 
stability characteristics of the jet are masked by addition of external forcing. Note that 
these simulations are of a transitional flow, whereas the experimental data are obtained in 
a fully turbulent jet flame. The imposition of explicit forcing was therefore necessary to 
mimic a turbulent environment so that meaningful comparisons with experimental results 
could be made. This forcing, however, overcomes the mixing suppression caused by factors 
such as reaction exothermicity, compressibility, etc. [14]. In order to make this point clear, 
simulations are performed of a heat releasing flame (Flame D). This flame has the same 
characteristics ( i.e . identical velocities, Damkolher, and Zeldovich numbers) as Flame A, 
but it also accounts for the effects of exothermicity ( e.g . a non-zero enthalpy of combus- 
tion). A comparison between the total products generated by this flame and that by Flame 
A is made in Fig. 9. This figure indicates a higher rate of product formation for the case 
of exothermic reaction. This observation is not consistent with any of the previous DNS 
[8, 14] or experimental [32] results. The reason for this discrepancy is due to the effects of 
the Arrhenius kinetics model in the present forced simulations, and cannot be observed in 
constant rate kinetics simulations and laboratory experiments [8, 14, 32]. The increased level 
of heat release also delays the onset of extinction, since the heat transfer away from the flame 
is countered by the heat generated inside. This new observation, therefore, suggests that in 
the setting of a “turbulent” flame in which the “transitional effects” are not dominant, and 
where the chemical reaction is describable by an Arrhenius model, the heat release results 
in a higher reactant conversion, even though the rate of mixing may be somewhat less. 
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Based on our observation, we recommend that the effects of exothermicity be assessed by 
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means of laboratory measurements. These measurements must involve a reacting system 
whereby the rate of reaction conversion is temperature dependent (unlike that employed in 
[32]), and in which the large scale mixing intensity is not significantly affected by the heat 
release (see [11, 12, 17, 24, 19, 18] for discussions on the effects of heat release in low speed 
combustion, and [14, 20] for discussions of such effects in high speed combustion). It is 
important to note that the issue would not be resolved by a mere comparison of the mixing 
characteristics between a reacting and a non-reacting layer. An appropriate analysis requires 
the consideration of two reacting systems with the same hydrodynamic characteristics but 
with different magnitudes of the enthalpy of combustion. The extent of validity of our 
conclusion under more complex chemistry models can be determined by these experiments. 


4 Concluding Remarks 

A two-four compact differencing algorithm is employed for direct numerical simulations of 
a spatially developing planar jet flame under non-equilibrium chemical conditions. With 
the high numerical accuracy provided by this algorithm, it is possible to perform direct 
simulations under physically realistic conditions of variable density and exothermic reaction. 
Here, the kinetic mechanism is assumed to obey the idealized model of the type F + 0 — ► 
Product + Heat. 

The data produced by DNS are statistically sampled to assess the compositional structure 
of the flame. The results of this analysis are shown to be consistent with those of laboratory 


16 


measurements, even with the assumption of the idealized chemistry model. It is shown that 
as the intensity of mixing increases, and the flame approaches extinction, the magnitudes 
of the ensemble mean and rms of the product mass fraction decrease whereas those of the 
reactants increase. This conclusion is only valid if the physical mechanisms by which the 
growth of instability waves is suppressed are masked. This procedure was implemented here 
by means of imposing a low amplitude random forcing at the entrance of the jet. 

The marginal and the joint PDF’s of the reactants mass fractions, extracted from DNS 
data, show features in accord with laboratory data. That is, as the relative intensity of 
mixing is increased, and the value of the Damkohler number is decreased, the PDF’s show 
a lesser probability of product formation. The scatter plots of the instantaneous product 
formation also show trends in accord with experimental data. These results also scale well 
with those based on the laminar diffusion flamelet model. However, the strong scatter in 
the results at the low values of the Damkohler number considered, suggest that to achieve 
a better agreement requires the use of a “flow dependent” form for the scalar dissipation. 
It is anticipated that the “correct” functional form of this dissipation is strongly dependent 
on the flow configuration, and it is not clear what form is appropriate for the setting of the 
unsteady jet flame considered here. 

The results further indicate that in the setting of a “turbulent” flame, the exothermicity 
results in an enhanced reactant conversion and an increased rate of product formation. This 
trend is not in accordance with previous DNS and laboratory results. The factors resulting 
in this discrepancy are: (1) imposition of external forcing, (2) the implementation of an 
Arrhenius kinetic model. Both of these factors are very important in diffusion flames in 
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which the flow is fully turbulent and the chemical reaction rate can be assumed to obey 
the Arrhenius law. Therefore, we propose that in typical hydrocarbon turbulent diffusion 
flames, the effect of heat release is to increase the rate of product formation even though the 
rate of mixing may be somewhat reduced. 

The results of our current preliminary investigations indicate that our major conclusions 
remain unchanged for 3-D simulations. However, the effects of small scale turbulence motion 
are very important in the analysis of the compositional structure with a view on flamelet 
modeling [26, 5, 6]. With such simulations, it is anticipated that the scatter would become 
more severe. Also, the extraction of high level statistical data, such as the marginal and 
joint PDF’s are more appropriate in the random field of an unsteady three-dimensional 
configuration, dominated by small scale structures. 

Considering the fact the kinetics mechanism employed in the present simulations is very ide- 
alized, the qualitative agreement between the experimental measurements and the simulated 
data is noteworthy. 
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Symbols 


C. Concentration. 

c p . Specific heat of the mixture at constant pressure. 

D. The jet width at the inflow. 

Da. The Damkoher number: Da = 

E. Activation energy for the Arrhenius chemical reaction. 

F. Fuel. 

Kj. Pre-exponential factor of the reaction rate. 

O. Oxidizer. 

Me. The convective Mach number. 

P. Product. 

V. Single-Point Probability Density Function. 

Pr. The Prandtl number. 

Q. Non-dimensionalized heat release parameter Q = 

It. Universal gas constant. 

Re. The Reynolds number. Re = ^ U P — . 

Sc. The Schmidt number. 

T. Temperature. 

U. Streamwise velocity. 

x,y. Streamwise and cross stream directions. 

Y. Mass fraction. 

Ze. The Zeldovich number, Ze = ^ . 
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Greek Symbols, Subscripts and Other Notations 


v. Kinematic viscosity. 
p. Density. 

( ). Ensemble-averaged. 

—AH 0 . Enthalpy of combustion. 

AU. Velocity difference of the two streams = Uf — Uo- 
X . Instantaneous dissipation of the conserved mixture fraction. 
ip. Representation of a reacting scalar variable. 
u. Reaction rate. 

uty. Appropriate reaction rate for the scalar variable rft. 

£. The mixture fraction. 

Tp. Total product mass fraction. Tp(x ) — ££ Yp{x,y)dy. 
st. Stagnation. 
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6 Appendix II 


This work summarizes the first efforts in validating the use of the Pearson family of frequen- 
cies in the modeling of non-homogeneous reacting flows. This work was deemed necessary 
for the following reasons: 

1. Before applying the Pearson density for subgrid modeling in LES, this closure must be 
validated and assessed for RANS of non-homogeneous flows. 

2. Because of the similarity of the PDF’s produced by the AMC and the Pearson family 
in homogeneous flows (Appendix I), it is speculated that the same similarity holds in non- 
homogeneous flows. However, the implementation of the AMC for non-homogeneous flow is 
very difficult. Therefore, in this work, we have started with the simple model of the Beta 
density. 

Work is in progress for the extension of this work for the LES of non-homogeneous shear 
flows. 
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7 Appendix III 


In this work, Mr. Craig Steinberger, with the help of his colleague, Mr. Thomas Vidoni, 
has focused his efforts in further investigations of non-equilibrium compressible reacting 
shear flows. This particular aspect of his work is concerned with: (1) The effects of finite- 
rated kinetics in statistical description of non-equilibrium diffusion flames, (2) The extent of 
validity and applicability of the laminar diffusion flamelet modeling for finite rate chemistry 
diffusion flames, and (3) Further effects of exothermicity in turbulent diffusion flames. 

In addition to the subject matters described above, the data bank provided by these sim- 
ulations would be very useful for further statistical modeling of turbulent flames in future 
investigations. 
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Figure Captions 
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Figure 1. (a) Schematic diagram of the planar jet configuration, (b) the grid. 

Figure 2 (a). Plot of product mass fraction contours, (a) Instantaneous data at several times 
(f 3 > t 2 > 1 1 ). (b) Ensemble- averaged data. 

Figure 3. The cross-stream variations of the ensemble-averaged product mass fraction at 
x — 10D. 

Figure 4. The cross-stream variations of the rms of the product mass fraction at x = 10 D. 

Figure 5. Joint PDF’s of the reactants mass fractions, ’P(Yf,Yo) at x = 10 D,y = 3.8 D. (a) 
Flame A. (b) Flame B. 

Figure 6. Marginal PDF’s of product mass fraction V(Yp) at x = 10 D,y = 3.8D. (a) Flame 
A, (b) Flame B. 

Figure 7. Scatter plots of product mass fraction vs. the mixture fractions for Flame A, and 
comparison with the Laminar Diffusion Flamelet Model, (a) x = 6.8 D, (b) x = 7.9 D, (c) 
x = 10D, (d) x = 11. 4D. 

Figure 8. Scatter plots of product mass fraction vs. the mixture fractions for Flame B, and 
comparison with the Laminar Diffusion Flamelet Model, (a) x = 6.8D, (b) x = 7.9 D, (c) 
x = 10 D, (d) x = 11.4D. 

Figure 9. The streamwise variation of the total product mass fraction 7p(x) under both 
non-heat releasing (Flame A) and heat releasing (Flame D) conditions. 
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Figure 7d 

Da = 0.50, i =210 
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Figure 8a 

Da = 0.25, i = 125 





Figure 8b 

Da = 0.25, i = 145 









Figure 8d 

Da = 0.25, i = 210 








